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We  study  the  equation  of  state  of  sodium  using  the  molecular  dy- 
namics  technique  whereby  the  classical  motion  of  a  system  of  ions  is 
solved  with  the  aid  of  computers.  The  interaction  potential  between  p 
pairs  of  sodium  ions  consists  of  coulomb  and  Born-Mayer  repulsion  terms 
and  an  effective  ion-ion  Interaction  derived  from  pseudopotential 
theory.  This  theory  includes  the  effects  of  electron  gas  screening, 
exchange,  and  correlation.  We  use  a  model  pseudopotential  with  param¬ 
eters  fit  to  experimental  low- temperature  data.  By  using  this  technique, 
we  are  able  to  begin  with  an  atomic  description  of  a  simple  metal  and 
proceed  to  calculate  its  macroscopic  thermodynamic  properties. 

We  calculate  equation-of -state  points  consisting  of  the  total  in¬ 
ternal  energy  in  volume  and  temperature  space.  For  our  study,  the 
volume  ranges  from  10Z  expansion  to  10Z  compression  of  the  normal  den¬ 
sity  and  the  temperature  ranges  from  0  to  600  Kelvin.  We  are  able  to 
calculate  directly  values  of  the  function  that  contains  the  enharmonic 
contributions  to  the  energy.  We  report  the  results  for  calculations  of 
solid  sodium  in  the  hexagonal  close-packed  (hep)  and  body-centered 
cubic  (bcc)  phases,  and  of  liquid  sodium. 

At  high  temperatures  the  molecular  dynamics  system  melts.  We  cool 
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the  liquid  sodium  back  to  low  temperatures  and  it  forms  a  metastable 
glassy  state  for  which  we  are  able  to  calculate  equatlon-of -state  points. 
We  study  the  dynamics  of  the  melt  transition  and  define  a  region  where 
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partial  melting  occurs.  The  upper  limit  that  we  place  on  the  melting 
temperature  is  consistent  with  the  observed  value  and  the  calculated 
heat  of  fusion,  diffusion  coefficient,  and  atomic  distributions  agree 
veil  with  experiment. 

We  illustrate  the  unique  capabilities  of  the  molecular  dynamics 
technique  by  inducing  a  dynamic  bcc-to-hcp  martensitic  phase  change. 

We  change  the  shape  of  the  calculations!  volume,  which  pushes  the  bcc 
sodium  structure  over  a  potential  hill.  It  then  spontaneously  trans¬ 
forms  to  the  more  stable  hep  structure. 

The  results  of  this  study  demonstrate  that  the  molecular  dynamics 
technique,  coupled  with  an  interaction  potential  that  adequately  de¬ 
scribes  the  ion- ion  interaction  In 
late  the  macroscopic  properties  of 
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We  study  the  equation  of  state  of  sodium  using  the  molecular  dy¬ 
namics  technique  whereby  the  classical  motion  of  a  system  of  ions  is 
solved  with  the  aid  of  computers*  The  interaction  potential  between 
pairs  of  sodium  ions  consists  of  coulomb  and  Born-Mayer  repulsion  terms 
and  an  effective  ion-ion  interaction  derived  from  pseudopotential 
theory*  This  theory  includes  the  effects  of  electron  gas  screening, 
exchange,  and  correlation.  We  use  a  model  pseudopotential  with  param¬ 
eters  fit  to  experimental  low- temperature  data*  By  using  this  technique, 
we  are  able  to  begin  with  an  atomic  description  of  a  simple  metal  and 
proceed  to  calculate  its  macroscopic  thermodynamic  properties. 

We  calculate  equation-of -state  points  consisting  of  the  total  in¬ 
ternal  energy  in  volume  and  temperature  space*  For  our  study,  the 
volume  ranges  from  10Z  expansion  to  10Z  compression  of  the  normal  den¬ 
sity  and  the  temperature  ranges  from  0  to  600  Kelvin*  We  are  able  to 
calculate  directly  values  of  the  function  that  contains  the  enharmonic 
contributions  to  the  energy*  We  report  the  results  for  calculations  of 
solid  sodium  in  the  hexagonal  close-packed  (hep)  and  body-centered 
cubic  (bcc)  phases,  and  of  liquid  sodium* 

At  high  temperatures  the  molecular  dynamics  system  melts.  We  cool 
the  liquid  sodium  back  to  low  temperatures  and  it  forms  a  metastable 
glassy  state  for  which  we  are  able  to  calculate  equatlon-of-state  points. 
We  study  the  dynamics  of  the  melt  transition  and  define  a  region  where 


partial  melting  occurs.  The  upper  limit  that  we  place  on  the  melting 
temperature  Is  consistent  with  the  observed  value  and  the  calculated 
heat  of  fusion,  diffusion  coefficient,  and  atomic  distributions  agree 
well  with  experiment. 

Ve  illustrate  the  unique  capabilities  of  the  molecular  dynamics 
technique  by  inducing  a  dynamic  bcc-to-hcp  martensitic  phase  change. 

We  change  the  shape  of  the  calculational  volume,  which  pushes  the  bcc 
sodium  structure  over  a  potential  hill.  It  then  spontaneously  trans¬ 
forms  to  the  more  stable  hep  structure. 

The  results  of  this  study  demonstrate  that  the  molecular  dynamics 
technique,  coupled  with  an  interaction  potential  that  adequately  de*» 
scribes  the  ion- ion  interaction  in  a  simple  metal,  can  be  used  to  calcu¬ 
late  the  macroscopic  properties  of  such  systems. 
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We  study  the  equation  of  state  of  sodium  using  the  molecular  dy¬ 
namics  technique  whereby  the  classical  motion  of  a  system  of  icns  is 
solved  with  the  aid  of  computers.  The  interaction  potential  between 
pairs  of  sodium  ions  consists  of  coulomb  and  Born-Mayer  repulsion  terms 
and  an  effective  ion-ion  interaction  derived  from  pseudopotential 
theory.  This  theory  includes  the  effects  of  electron  gas  screening, 
exchange,  and  correlation.  We  use  a  model  pseudopotential  with  param¬ 
eters  fit  to  experimental  low- temperature  data.  By  using  this  technique, 
we  are  able  to  begin  with  an  atomic  description  of  a  simple  metal  and 
proceed  to  calculate  its  macroscopic  thermodynamic  properties. 

We  calculate  equation-of -state  points  consisting  of  the  total  in¬ 
ternal  energy  in  volume  and  temperature  space.  For  our  study,  the 
volume  ranges  from  10%  expansion  to  10%  compression  of  the  normal  den¬ 
sity  and  the  temperature  ranges  from  0  to  600  Kelvin.  We  are  able  to 
calculate  directly  values  of  the  function  that  contains  the  anharmonic 
contributions  to  the  energy.  We  report  the  results  for  calculations  of 
solid  sodium  in  the  hexagonal  close-packed  (hep)  and  body-centered 
cubic  (bcc)  phases,  and  of  liquid  sodium. 

At  high  temperatures  the  molecular  dynamics  system  melts.  We  cool 
the  liquid  sodium  back  to  low  temperatures  and  it  forms  a  metastable 
glassy  state  for  which  we  are  able  to  calculate  equation-of -state  points. 
We  study  the  dynamics  of  the  melt  transition  and  define  a  region  where 
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partial  melting  occurs.  The  upper  limit  that  we  place  on  the  melting 
temperature  is  consistent  with  the  observed  value  and  the  calculated 
heat  of  fusion,  diffusion  coefficient,  and  atomic  distributions  agree 
well  with  experiment. 

We  illustrate  the  unique  capabilities  of  the  molecular  dynamics 
technique  by  inducing  a  dynamic  bcc-to-hcp  martensitic  phase  change. 

We  change  the  shape  of  the  calculational  volume,  which  pushes  the  bcc 
sodium  structure  over  a  potential  hill.  It  then  spontaneously  trans¬ 
forms  to  the  more  stable  hep  structure. 

The  results  of  this  study  demonstrate  that  the  molecular  dynamics 
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I.  INTRODUCTION 

Psfcwi© potential  theory  has  proved  useful  in  studying  many  prop¬ 
erties  of  simple  metals  (for  example,  see  Refs.  1-7).  It  is  a  method 
for  solving  the  Schrodinger  equation  which  contains  the  essential  fea¬ 
tures  of  the  behavior  of  the  electrons  in  these  metals. 

Molecular  dynamics  is  a  technique  for  studying  the  classical  be¬ 
havior  of  a  many-particle  system.  Newton’s  second  law  is  solved  from 
the  force  between  pairs  of  interacting  particles  where  the  force  is  de¬ 
termined  by  the  gradient  of  the  interaction  potential  between  pairs  of 
particles.  The  interaction  potentials  most  commonly  used  in  molecular 
dynamics  calculations  are  empirically  determined  (for  example,  the 
Lennard-Jones  potential) .  While  such  potentials  are  easy  to  use  and 
adequately  represent  the  behavior  of  some  systems,  they  are  not  appro¬ 
priate  to  the  interactions  between  the  ions  in  a  simple  metal. 

We  propose  to  use  an  effective  ion- ion  interaction  potential  de¬ 
rived  using  the  pseudopotential  method  in  our  molecular  dynamics  calcu¬ 
lations.  By  doing  this  we  are  able  to  start  with  the  atomic  description 
of  the  simple  metal  and  proceed  to  calculate  macroscopic  thermodynamic 
properties.  This  effective  interaction  potential  is  long  range  and  re¬ 
quires  the  inclusion  of  approximately  170  neighbors  for  each  particle 
when  the  forces  are  calculated. 

Additionally,  the  crystal  potential  energy  of  the  effective  inter¬ 
action  is  only  a  fraction  (about  (0.2Z)  of  the  total  energy  of  the 
crystal.  The  volume-dependent  energy  terms  (such  as  the  electron-gas 
kinetic  energies  and  the  exchange  and  correlation  energies)  are  prima¬ 
rily  responsible  for  holding  the  crystal  together.  These  terms  are 
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calculated  separately  and  added  to  the  structure-dependent  energy,  which 
is  calculated  by  the  molecular  dynamics  program. 

In  the  "Background"  section  of  this  paper  we  describe  the  theory 
for  both  the  pseudo potential  method  and  the  molecular  dynamics  technique. 
This  discussion  has  been  developed  elsewhere,  as  referenced  in  the  text, 
and  we  include  the  features  important  to  our  study  for  completeness. 

In  Sec.  Ill  and  IV  we  describe  the  calculation  of  the  effective 
ion- ion  interaction  for  sodium  and  the  setup  of  our  molecular  dynamics 
calculations.  We  outline  in  detail  the  necessary  steps  for  determining 
the  parameters  and  run  conditions  necessary  to  perform  these  calculations. 

In  Sec.  V  we  develop  the  equations  necessary  for  the  calculation  of 
the  volume-dependent  energy  terms,  mentioned  above,  that  must  be  added 
to  the  structure-dependent  terms  to  arrive  at  the  total  system  energy. 

In  these  first  five  sections,  we  describe  the  technique  for  using 
the  quantum  mechanical  results  of  pseudopotential  theory  in  a  classical 
trajectory  calculation  of  the  motion  of  the  ions  of  sodium.  Using  this 
technique,  we  calculate  equation-of-state  points  consisting  of  the  sys¬ 
tem  energies  at  given  volumes  and  temperatures  for  solid  sodium  in  the  hep 
and  bcc  phases  and  for  liquid  sodium.  We  are  able  to  define  the  equa¬ 
tion  of  state  for  a  "glassy"  sodium,  which  is  the  liquid  extended  to  a 
metastable  low- temperature  state.  The  results  of  these  calculations 
are  direct  measurements  of  the  anharmonic  contributions  to  the  system 
energy  for  these  systems. 

One  of  the  major  advantages  of  the  molecular  dynamics  technique  is 
the  ability  to  follow  the  dynamics  of  the  system  studied.  We  discuss 
the  melting  of  sodium  and  are  able  to  define  a  transition  region  of 
partial  melting.  The  calculated  heats  of  fusion  and  upper  limit  to  the 


3 

melting  temperature  are  in  agreement  with  observed  values.  We  are  also 
able  to  calculate  the  atomic  distributions  and  diffusion  coefficients 
for  sodium. 

As  a  final  result  and  as  a  demonstration  of  the  unique  capabil¬ 
ities  of  molecular  dynamics  calculations,  we  present  a  calculation  of 
the  dynamics  of  the  bcc->*hcp  phase  change  in  sodium.  We  artificially 
change  the  shape  of  the  calculational  volume.  This  pushes  the  meta- 
stable  low-temperature  bcc  structure  over  a  potential  hill.  The  system 
is  then  able  to  spontaneously  relax  into  the  preferred  hep  structure 
with  an  accompanying  increase  in  temperature. 

Our  results  indicate  that  the  technique  described  in  this  paper 
adequately  represents  the  behavior  of  sodium  in  the  volume  (10%  expan¬ 
sion  to  10%  compression)  and  temperature  (0-600  K)  ranges  studied. 
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II.  BACKGROUND 

This  section  discusses  the  necessary  theoretical  background  for 
development  of  the  technique  described  in  the  introduction,  by  which 
we  will  study  Che  sccium  equation  of  state.  Section  II. A  discusses  the 
pseudo potential  method  and  culminates  with  an  expression  for  the  total 
effective  ie*ri-ion  interaction.  Section  II. B  describes  the  molecular 
dynamics  tech*.  and  the  applicable  equations. 

A.  The  Pseudopotential  Method 

Wr  calculate  the  ion  motion  in  a  simple  metal.  A  simple  metal  is 

one  for  which  the  conduction  electrons  behave  very  nearly  as  if  they 

1  2 

comprise  a  free-electron  gas.  ’  We  assume  that  the  metal  consists  of 

ions  and  conduction  electrons.  The  ions  are  composed  of  the  nuclei  and 

the  core  electrons.  An  ion  core  does  not  overlap  with  other  ion  cores. 

The  core  electron  states  are  assumed  to  be  the  same  as  the  respective 
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states  in  a  free  atom.  9 

In  this  section  we  describe  the  pseudopotential  method,  which  is  a 
technique  for  solving  the  Schrddinger  equation  for  the  energy  of  the 
conduction  electrons.  In  Sec.  II. A. 1  we  express  the  electron  energy  in 
terms  of  the  pseudopotential.  In  II. A. 2  we  restructure  the  equation 
to  express  the  electron  energy  in  terms  of  an  effective  ion-ion  inter¬ 
action.  In  Secs.  II. A. 3  -  II. A. 5  we  incorporate  a  local  approximation  to 
the  pseudopotential  and  modify  the  theory  to  include  electron  screening, 
exchange,  and  correlation  effects.  In  Sec.  II. A. 6  we  develop  the  model 
pseudopotential  that  we  will  use  for  this  study,  and  in  Sec.  II. A. 7  we 
add  the  coulomb  and  core  repulsion  terms  to  obtain  the  total  ion-ion  inter¬ 
action,  which  will  be  used  in  the  molecular  dynamics  calculations. 
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1.  Conduction  Electron  Energy.  The  total  Hamiltonian,  for 

a  system  of  N  ions  and  N’  electrons,  neglecting  external  interactions, 
is 

“tot  ■  al  *  a.  <» 


with  subscript  e  denoting  electrons  and  subscript  I  denoting  ions, 
where 


“i  »!<*»- *k> 

l  1  l  k>£, 


T  +  V 
I  I 


1 } 


■z  ±  +ES  v?i  -  'V  +ZZ 

i  e  i  j  >i  i  l 


T  +  V  +  VT 
e  e  I-e 


£,k  -  1,  •  •  *N 
i,j  -  1  * • • *NT 
r^  ■  electron  positions 
R^  *  ion  positions 
T  *  kinetic  energy 
m  *  mass 

Vj  *  ion-ion  potential  interaction 
Ve  *  electron-electron  potential  interaction 
■  ion-electron  potential  interaction 

The  Schrodlnger  equation  is 


-  E 


w 

TOT  “TOT 


(2) 


i 


1 


^OT  TOT 
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The  solution  of  this  equation  is  intractable*  and  we  simplify  it 
with  another  well-known  assumption — the  adiabatic,  or  Bom-Oppenheimer , 
approximation.  The  essence  of  this  approximation  is  that  the  electrons 
readjust  themselves  so  rapidly  to  a  change  in  ion  configuration  that 
the  ions  are  regarded  as  fixed  when  solving  for  the  electron  energies. 
This  uncouples  the  electron  part  of  the  equation  from  the  ion  part. 
Therefore,  the  total  wave  function  is  separable  and  is  written  as 

^TOT  W  ' 


We  may  solve  the  electron  problem  for  fixed  ion  positions  to  obtain 


H  y 
e  e 


E  V 
e  e 


(3) 


This  substituted  in  Eq.  (2)  yields 


Wtot  ■  «t  *  VI  +  We  ■  W, 


(4) 


where 


TJi'  ¥  +  VJFJF  +  E 

lie  lie  ele 


EV  „  h2  a2 

JL.  v  f  .  V  - - 2 -  f  if 

2m  TI  e  2m  ~  T~T 


*• 


— .  *.2  s2  h2  __  /  a?  3yT  a2'!'  \ 


l  l 


The  last  two  terms  in  this  expression  are  neglected  in  the  adiabatic 

approximation  because  they  contribute  negligibly  to  the  system  en- 
A  5 

ergy.  *  Thus,  we  can  write 


3 

■«*? 


v4 


or 


Sm'tor  ■  VTI  *  VI  +  VTI  ■  WA 


> 'v  i  *  ■ 
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(II,  +  ».)»,  -  EmTy,  . 


and  we  see  that  the  adiabatic  approximation  allows  consideration  of  the 

electron  problem  separately  from  the  ion  problem.  The  conduction  elec- 

3 

tron  energy  becomes  an  effective  potential  energy  of  the  ions.  We  now 
solve  the  electron  problem. 


HF  «Ef 
e  e  e  e 


We  assume  the  self-consistent  field,  or  one-electron,  approximation, 
where  the  potential  (V(r))  that  an  electron  moves  in  is  calculated  in  a 
self-consistent  manner.^  We  write 


H  -  T  +  V(r)  , 
e  e 


where  V(r)  is  the  self-consistent  field  and  we  write  the  Schrodinger 
equation  for  a  single  electron  as:^ 


He^  *  fTe  +  *  £e^ 


Following  the  notation  in  Ref.  6,  we  use  the  index  t  to  denote  core 
electron  states  and  j  to  denote  the  state  centered  at  ion  position  R j . 

We  have 

[T,  ♦  .  (6) 

where  the  ^  .  are  the  same  as  for  the  free  atom,  according  to  our  first 

t * j 

assumption. 

With  the  above  equation  established  within  the  constraints  of  the 
assumptions  mentioned,  we  now  restructure  it  in  terms  of  the  pseudo¬ 
potential  formulation.  First  we  expand  the  electron  wave  functions 
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in  terms  of  a  basis  of  functions  that  are  constructed  to  be  orthogonal 
to  the  core  states.  Following  Harrison,1’*’  we  use  the  notation 


and 


h  ik-r 


\k>  =  Q  ^  e 
I  t,j>  =  ’■l>t (r  -  R  ) 


,  a  normalized  plane  wave, 

,  a  normalized  core  function 
centered  at  ion  position  Rj , 


<t,j |k> 


/ 


d3r 


^t(r 


*  v  ik-r 

Rj)e 


Then  we  write  the  basis  functions,  called  orthogonal ized  plane  waves,  as 


opwk  -  (i  -  p) |£>  , 


where  P  is  the  projection  operator,  which  projects  any  function  onto  the 
core  states1  (note  that  <t,jjtf,jf>  -  <5^,5^,), 

p  |t,jxtj|  .  (7) 

j 


We  expand  the  conduction  electron  wave  function  in  terms  of  the  orthog- 
onalized  plane  waves  to  obtain 

*  -  (1  -  P)]T  ak  1  k>  ,  (8) 

k 

where  a^  are  the  expansion  coefficients.  If  we  substitute  in  the 
Schrodinger  equation,  the  solution  could  be  attempted  using  the  standard 
orthogonalized  plane  wave  method.  We  are  attempting  a  different  approach 
and  will  restructure  the  equation.  We  introduce  the  pseudowavefunction, 

“  53  a^|k>  (9) 

so  that 


(10) 
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v. 


We  note  that  $  is  an  expansion  of  free  electron  states  and  that  <p  ■  ^ 
outside  the  cores  because  the  projection  operator  is  zero  there. ^  Sub¬ 
stituting  in  the  Schrodinger  equation  yields 

Te(l  -  P)0  +  V(r)  (1  -  P)<j>  -  E(1  -  P)<j>  , 

and  rearranging  gives 

Te4>  +  V(r)<j>  -  [Te  +  V(r)]P$  +  EP<t>  «  E$  , 

so  that 


T  <p  +  =  E$  ,  (11) 

e 

where  we  have  defined  the  operator, 

W  *  V(r)  -  [T^  +  V(r)]P  +  EP  =  the  pseudo potential  operator. 

Using 

[Te  +  V(r)]P  EtJ|t,jxt,j|  , 

we  have 

w  -  V(r)  +  £  (E  -  Et>J)|t,jxt,j|  ,  (12) 

t,j 

which  is  the  pseudopotential  equation.^ 

As  Harrison1’^  notes,  the  pseudo potential  has  several  interesting 
properties.  It  is  nonlocal  in  that  it  depends  on  all  ion  positions  and 
states.  It  is  an  operator  and  is  not  restricted  to  multiplying  the 
wave  function.  Its  form  is  not  unique  and  the  pseudowavef unctions  are 
not  unique.  That  is,  an  arbitrary  number  of  core  states  can  be  added 
to  the  pseudowavef unctions  and  they  will  still  solve  the  Schrtidinger 
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equation.  Also,  Harrison  shows  that  (E  -  E  )  may  be  replaced  by  any 

t  f  j  ' 

function  of  the  energies  and  the  core  states  [f(E,t,j)]  and  the  solution 

will  remain  unchanged.  There  exist  many  valid  forms  that  will  yield 

correct  energies  and  wave  functions.** 

The  pseudopotential  property  of  Importance  here  is  that  it  can  be 

considered  small.  V(r)  is  negative,  P  is  positive,  and  (E  -  E  )  is 

t » J 

positive  so  that  the  terms  in  Eq, (12) tend  to  cancel.^  Experimental 
evidence  corroborates  that,  in  simple  metals,  the  conduction  electrons 
behave  much  like  free  electrons. 

Given  that  W  is  small,  we  may  use  perturbation  theory  to  solve  for 
the  electron  energies.  We  note  here  that,  even  though  the  W  form  is 
arbitrary,  if  the  problem  were  solved  exactly  (i.e.  jo  all  on  in  W)  , 
the  correct  solution  would  always  be  obtained.  However,  using  perturba¬ 
tion  theory,  the  W  form  will  affect  the  result  and  the  validity  of  the 
form  used  must  be  determined  by  the  results  obtained  when  applied  to  a 
specific  problem.^ 

Using  perturbation  theory,  we  calculate  the  electron  energy  in  the 
state  k  to  the  second  order, ^ 


\  -  Ek  +  <siw|£>  +  y; 


* <k|w|k  +  q><k  +  q|wlk> 


£k  "  £k-h, 


(13) 


where 


nV 

2m 


(U) 


are  the  free-electron  kinetic  energies.  The  prime  on  the  summation  in¬ 
dicates  that  the  q-0  term  is  omitted  from  the  sum.  We  also  may  calcu¬ 
late  the  first-order  pseudowavefunctions  as 


i : 


■  |ts>  +  Saqklk  +  q>  * 


where 
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(15) 


(16) 


The  q  -  0  term,  which  does  not  Immediately  concern  us,  is  determined  in 
Ref.  2  by  normalization. 

Given  an  appropriate  pseudo potential,  we  evaluate  this  expression 
and  sum  over  all  occupied  electron  states  to  determine  the  electron 
energy. 


(17) 


nk  *{o 


It  is  occupied 
k  is  not  occupied 


This  particular  result  helps  solve  the  ion  motion  problem  (see 
Eq.  5).  In  evaluating  the  conduction  electron  energy  we  must  consider 
the  electron  self -energy  terms  that  are  counted  twice  in  Eq.  (17) .  This 
will  be  discussed  appropriately  in  the  ensuing  sections. 

2.  Effective  Ion- Ion  Interaction.  We  evaluate  now  the  matrix  ele¬ 
ments  in  the  perturbation  expansion  for  the  electron  energies  [Eq.  (13)]. 
In  Eq.  (12)  and  its  discussion,  the  general  expression  for  the  pseudo¬ 
potential  was  shown  to  be 


w  -  V(?)  +  £  f(E,t,j)|t,jXt,j|  . 


(18) 


V(r)  contains  the  potential  field  due  to  the  ions,  V^(r),  which  may  be 
written  as  the  sum  of  contributions  from  individual  ions1  at  positions 

rr 


tej/ 31**.%:* .  t  '■  > '  *  v  '•  w » n* n  &  * 
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;i&  ■£  vi(? "  fy  • 


The  ion  screening,  according  to  Harrison  may  be  superimposed  at  the 
ion  sites  and  will  be  spherically  symmetric.  We  have  already  written 
the  core  states  as 


|t,j>  -  i|>t(r  -  . 


If  we  require  that  the  function  f(E,t,j)  depends  only  on  the  core  states 
through  the  index  t,  we  may  separate  the  pseudopotential  and  write  it 
as  a  sum  of  contributions  from  individual  ion  sites  as^ 


w  -y;  w(?  -  . 


This  assumption  is  essential  to  the  pseudopotential  matrix  element  cal¬ 
culation. 

Now  we  are  able  to  factor  the  matrix  elements.  Choosing  W  of  Eq. 
(20)  as  our  pseudopotential  allows  factoring  out  the  structure-dependent 
term.  We  write 


+  ,  1  C  -i(£+q)  rr  x  i«.?  .3 

<k  +  q|W|k>  *  —  J  e  M  w(r  “  rj)e  ^  " 


-iq.rj 


Changing  Che  summation  and  integration  order  and  factoring  out  e 


<$  +  q  |  W  |  k>  -  e~iq*r  j  J  e"i(&3)  •  (?-?j)w(?  .  J 


The  integral  is  just  the  integration  of  an  individual  potential  with 
respect  to  the  position  of  that  particular  ion.  The  sum  contains  the 
information  about  the  system  structure,  so  we  write 


wt*  -a4-!*-- 


(21) 


S (q)  -  structure  factor  =  —  elq  J  , 

i 


and  we  define 


<ic  +  qjw|k>  *  form  factor  =  r  w(r)e^*^ 


(22) 


where 


■  atomic  volume  ■  &/N 


So  the  pseudopotential  matrix  element  becomes 


<k  4*  q|w|k>  -  S(q)<k  +  qjw|k>  . 


(23) 


This  factoring  is  critical  to  the  development  of  the  pseudo potential 
method.  It  allows  simple  solution  of  the  seemingly  intractable  many- 
body  problem.1  The  detailed  ion  positions  enter  only  through  the  struc- 
ture  factor  S(q)f  and  the  ionic  potential  details  enter  only  through  the 
form  factor. 

We  are  ready  to  sum  electron  energies  over  the  available  states  to 
determine  the  total  electron  energy.  We  rewrite  Eq.  (13)  In  terms  of  the 
structure  and  form  factors  as 

e.  ■  e.  +  s(o)<£|«|£>  +T'  tl&I&dkill  ±  +  q|»|t>  .  (24) 

R  q  ek  ek+q 

Note  here  that 

N 

s(o)  - IZ1 " 1  • 

j-i 

For  free  electrons  in  the  ground  state,  the  available  energy  states 


are  described  by  the  Fermi  sphere  of  radius 
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0  x/  j 

kf  "  fe" ) 


where  Is  the  Fermi  wavenumber.  With  the  pseudo  potential  present,  the 

energy  surface  is  not  spherical.  However,  when  evaluating  the  energies 

1  2 

to  second  order,  we  may  neglect  this  higher  order  effect  *  and  sum 
over  the  Fermi  sphere.  We  divide  by  N  to  obtain  the  total  electron 
energy  per  ion,  N 


iL 


t  -  r  . 

K  (2ir)Jy0 


where  we  have  used  the  density  of  states  in  wave-number  space  to  convert 
the  sum  to  an  integral,  including  a  factor  of  2  to  account  for  spin 


states . 


We  now  evaluate  the  contributions  to  Eq.  (26)  of  the  three  terms  in 


Eq.  (24).  The  first  is 


2n0  rkf  ,3.  2n0  /*k f  ft2k2  ,3, 

— 3j  ekd  k  m7ZJj  ^i"dk 

(2tt)  •'q  (2tt)  jq 


!jsk?-zK  • 


which  is  the  average  kinetic  energy  of  the  electrons  times  the  valence. 
The  second  term  is 


2n0  /*5: 
(2  it  )3Jq 


f  f<k|v|k*l 
Jq 


2!i0  kf:/^£HfcA 

^  3  /  k  d  k 

«*>  •£  f0  'ih 

.  3-  V  4  •  2  <?H*>  , 
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which  Is  2  times  the  average  value  of  <k/w/k>.  Both  of  these  terms  de¬ 
pend  on  volume  but  not  on  the  details  of  the  ion  positions,  and  they 
represent  the  free-electron  energy, ^  Exchanging  the  order  of  the  sum 
and  integral,  the  third  term  is 


Here,  we  define  the  part  of  this  expression  that  is  a  function  of  |q| 
only  as  the  energy  wave-number  characteristic  F(q).  It  is  determined 
by  w,  which  is  spherically  symmetric,  and  k^. 


F(q) 


(29) 


The  third  term  is  called  the  band-structure  energy,  E^g,  and  is  written 

Ehg  S*(?)S(?)F(q)  .  (30) 

q 

Note  that  F(q)  depends  on  the  volume  but  not  on  the  detailed  ion  arrange¬ 
ment,  which  is  determined  by  the  structure  factor.  This  band  structure 
energy  Interests  us  when  we  calculate  the  effects  of  an  ion  position 
change  in  a  constant  volume  situation.  It  may  be  considered  an  effective 
(indirect)  interaction  between  ions.^* 

To  show  this  in  a  more  direct  manner  we  restructure  the  expression 

,or  %.• 

•£'  S*<,~)s«><,)  -L  ?(,)££  .-If  trVfj)  . 

q  q  N  i-lj»l 
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Note  that 


I-j  i-J"1 


so  chat 


Eb,  -E’i  F(’>(H  +£ 

q  N  '  i.j 


where  the  prime  on  the  summation  indicates  exclusion  of  the  i  «  j  term. 
Therefore, 


h.  -aZ'  • 

i.j  q 


(31) 


The  second  term  here  is  volume  dependent  and  we  define 


w?)  -|£’F(q)el*‘?  •  (32) 

q 

We  have  written  the  structure-dependent  terms  as  a  summation  over  the 
ions  of  an  effective  potential  that  is  dependent  only  on  the  positions 
of  the  ions  through  (r  ^  -  r^) .  is  the  effective  interaction  between 

the  ions  for  which  we  have  been  looking  and  which  we  now  must  evaluate. 
Replacing  the  summation  of  Eq.  (32)  by  an  integral. 


E 


a_ 

(2ir)3 


we  evaluate  the  angular  parts  of  the  integral  over  q  space  to  obtain 
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V 


vim>(?)  ’  7jC*f(<1)  5VIL  "2  d’  •  <33> 

This  is  a  two-body,  central  force  interaction  between  ions  that  may  be 
added  to  the  direct  ion- interact ion  terms  when  calculating  the  ion  mo¬ 
tion. 

3.  The  Local  Approximation.  Anticipating  that  we  will  approxi¬ 
mate  the  pseudopotential  by  a  real  model  potential  with  adjustable 

parameters,  w(r),  we  make  the  "local  approximation,"  which  will  simplify 

2 

the  rest  of  the  development.  We  assume  that  w  is  Independent  of  q  and 
we  neglect  its  nonlocal,  operator  character.  Therefore*  we  can  write 

<£  +  qjw(r)  |k>  m  f  w(r)e  iq*r  d3r  -  w  .  (34) 

Uq  J  4 

Additionally,  the  true  wavefunction  now  equals  the  pseudowavefunction 

2 

to  first  order  within  a  normalization  constant  and  can  be  written 

-  lk>  +53’  aqlt|k  +  q>  .  (35) 

q 

We  now  evaluate  the  energy  wave-number  characteristic  in  the  local 
approximation.  w(q)  is  taken  out  of  the  integral  of  Eq.  (29)  so  that 


2 

The  Integral  is  evaluated  over  the  Fermi  sphere  as 


2,  2h2  2 

-tr  y  n 


me 


i 


(e(q)  -  1)  , 


(36) 
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where 


and 


e(q)  -  1 


me2  /l  -  n2 

2*kfh2n2'  211 


l  jin 

l  -n 


which  yields^ 


(37) 


-Ojq  2 

F(q)  -  — V|wnp[£(q)  -  1]  •  (38) 

8  ire2  q 

e(q)  is  the  static  Hartree  dielectric  function  for  free  electrons.  It 
has  a  logarithmic  singularity  at  q  ■  2k^,  which  will  affect  signifi¬ 
cantly  the  calculated  interaction  potential  form. 

Constructing  the  local  pseudopotential  by  following  the  notations 
of  Ref.  2,  we  write 


W  -  W-  +  W  +  W  ,  (39) 

B  s  x 

where  W_  is  the  "bare- ion"  pseudopotential,  the  local  potential  by  which 

D 

the  electrons  interact  with  ions.  Wg  is  the  Hartree  screening  contribu¬ 
tion,  which  includes  only  the  coulomb  interactions  with  the  other  con¬ 
duction  electrons  determined  self-consistently.^  contains  contribu¬ 
tions  due  to  exchange  and  correlation.  will  be  chosen  later  as  an 
appropriate  model  potential.  Wg  is  determined  accurately  in  the  next 
section.  W^  must  be  approximated  and  is  discussed  in  Sec.  II. A. 5. 

4.  Self-Consistent  Screening.  Within  the  adiabatic  approximation, 
we  treat  the  conduction  electrons  as  though  responding  in  an 


* 
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electrostatic  sense  to  the  total  potential  W.  The  self-consistent 

electron  contribution,  Wg,  due  to  the  resulting  charge  density,  d(r)  may 

be  calculated  using  Poisson’s  equation.  (We  are  following  Wallace’s 
2 

method  • ) 


^2$(r)  -  -A-rre  d(r)  , 


where 

Wg(r)  -  e$(r) 

and 


so  that 


d(r)  *  electron  density 


^Wg(r)  -  -4Tre2d(r) 


We  now  expand  Wg(r)  and  d(r)  in  Fourier  series 


(40) 


^2w  (r)  -  y  W  eiq  r  -  -  y  q2W  eiq‘r  -  -4ire2  d(r) 

q  q 

■  -47re2  7  d 

q  q 

The  subscript  q  denotes  the  qth  Fourier  coefficient.  Thus  the  expan¬ 
sion  coefficients  for  Wg  and  d  are  related  through 

"m  -  (fl • 

To  calculate  the  screening  contribution  we  must  calculate  the  charge 
density,  remembering  that  we  are  calculating  the  electron  energies  to 


AT 

* 

3 


I 

I 

I 

I 
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I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


20 


second  order  and  the  vavefunction  to  first  order  in  perturbation  theory. 
The  charge  density  is  given  by 


d(?)  ■  l  > 

k 


(42) 


where  the  sum  is  over  all  electron  states  and  n^  is  the  occupation  number, 

2 

which  is  one  when  the  state  is  occupied  and  zero  when  unoccupied.  Using 

Eq.  (35),  which  is  the  expression  for  ^  accurate  to  first  order,  we 
2 

calculate  ,  keeping  only  first-order  terms, 


►X  -  n'l|'  * 


(43) 


We  identify  the  Fourier  coefficient  d^,  defined  by 


d(«  - 1  d  .1’-r 

q  q 

as 


ft”1  l  n 


k  -  o_1nz  -  n”Lz 


2fi"  £  Vqk  ’  q  *  0 


(44) 


Substituting  the  value  for  a^,  yields 


mW 


d3k 


■’  ■  '  iv/  k2  -  |J  +  • 


(45) 


where  we  have  used  the  density  in  wave-number  space  with  a  factor  of  2 
for  spin  to  convert  the  sum  to  an  integral  and  the  zero  order  electron 
energy. 


k  “  2m  k 


(46) 


1  f ,  ■  r.u  1  *  Vfc,  >  t  >■&***  *i  MW&  pH# 
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We  can  calculate  the  Fermi  energy  E^  and  Fermi  wave  number  k^  to 

2 

second  order  due  to  the  local  pseudopotential  perturbation.  This 
exercise  results  in  zero-order  terms  for  and  k^,  followed  by  second- 
order  corrections.  First-order  corrections  do  not  appear.  This  result 
allows  us  to  take  the  integral  /  dJk  over  the  zero-order  Fermi  sphere 
and  ensures  that  the  results  will  still  be  good  to  first  order.  We 
evaluate  the  integral 

f  d3k 
k2  -  |k  +  q  |2 

2 

over  the  Fermi  sphere  as  we  did  in  Sec.  1.A.3  and  calculate 


q"W 

d - 1  [1  -  e(q)]  ,  (47) 

q  4ire2 


where  £(q)  is  the  static  Hartree  dielectric  function  given  by  Eq.  (37). 

Substituting  into  Eq.  (41),  we  find  the  result  for  the  Fourier  co¬ 
efficients  of  W  is 
s 


wo„  »  w  fl  -  e(q)j  . 
sq  q 


Momentarily,  we  are  neglecting  W^  of  Eq.  (39),  so 


w  *  wn  +  w  , 

q  Bq  sq 


and  we  see  that 


(48) 


*  V/e00  *  (**> 

q  Bq 

so  that  £(q)  plays  the  role  of  the  dielectric  function. 1 

Now  that  we  are  including  the  coulomb  energy  of  the  conduction 
electrons,  ve  must  reevaluate  the  expression  for  the  structure-dependent 
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part  of  the  total  band  structure  energyr  Eq.  (30).  In  this  equation  we 
have  double-counted  the  electron  self-energy  Ege  and  must  subtract  it 
to  get  the  correct  result.  That  is, 

S,  •  (5C 

q 

The  prime  on  F’(q)  indicates  evaluation  of  a  more  appropriate  energy 
wave-number  characteristic,  which  can  be  used  in  Eq.  (33)  for 

The  electrostatic  electron  self-energy  per  ion  is  given  by 

I  E*e  *  4/A  d(?)  V?)  '  (5: 

which  (neglecting  the  q  *  0  term)  can  be  written  in  terms  of  Fourier 
coefficients  as 


so  that 


I  E  »  —  f  d  W 

N  ee  2N  £  q  s-q 


Wsq“(q2)dq 


-  E  -  I  q  w  W 

N  ee  8^e2  q  Sq  S"q 


From  Eq.  (20),  we  have 


W(r)  -  l  w(|r  -  r  |)  , 

i 


and  we  can  show  that 


(52) 


(53) 


(54) 


W  -  S(q)w  , 

q  q 


(55) 


where 


Wq  “  Tj/W(r)e’iq'r  d3r  * 
Wq  *  IT  J w^r^e  lq  r  d3,r  * 


and ,  as  in  Eq .  ( 21 ) , 


s<’> 


which  yields 


T  s  S  q2  |w 
a  -a 


2  I  1 2 


N  ee  an  2  "  <1  -q  sq 


We  may  substitute  this  into  Eq.  (50),  using  the  relations  between  W^,  Wg^f 


and  W_  to  arrive  at 
Bq 


EL  -y*’  s  S  M-|W  ]2  £ia2  -  1 

Ss  ^  qs-q  giTe2l  Bq  e(q) 


The  appropriate  energy  wave-number  characteristic  for  use  in  Eq.  (33) 
to  calculate  is,  therefore. 


F(q)  .foL|w  |2  SisLz. 

(q)  Bq'  e(q) 


This  expression  accounts  for  the  self-consistent  screening  of  the  elec¬ 
tron  gas,  within  the  order  of  our  perturbation  calculation. 

5.  Exchange  and  Correlation.  We  include  approximately  the  exchange 
and  correlation  effects  in  the  screening  calculation.  We  follow  the 


V, 


discussion  in  Ref.  2.  The  screening  calculation  outlined  in  the  preced¬ 
ing  section  yields  the  appropriate  result  with  the  inclusion  of  an  addi¬ 
tional  interaction  due  to  exchange  and  correlation,  W^,  so  that 


W  -  VL  +  W  +  W 
q  Bq  sq  xq 


We  assume  that  there  exists  an  average  interaction  I(q)  within  the  elec- 
2 

tron  gas. 


W  -  I(q)  d 
xq  q 


Then,  using  Eq.  (41), 


W  +  W  -  +  I  (q)  ]  d 

sq  xq  l  q2  J  « 


[i  -  f(q)]  «*  * 

q  q 


where 


f(q)  *  Kq)  * 

47te^ 


and  we  parallel  the  screening  calculation  algebra  to  obtain 


q  1  +  [e(q)  -  1][1  -  f (q) ]  * 


We  recalculate  the  electron  self-energy  appropriately  [seeEq.  (51)]  from 


o 

E  -  ~  >  dq(W  +  W  ) 
ee  2N  s-q  x-q 


and  reevaluate  the  energy  wave-number  characteristic.  This  all  follows 
the  development  in  the  preceding  section,  with  the  additional  [1  -  f (q) ] 
term  to  keep  track  of.  The  result  is 
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F(q) 


2  e(q)  -  1 _ 

1  +  Ie(q)  -  1][1  -  f(q)l 


(65) 


Heine^  discusses  the  formalism  that  we  have  just  outlined  and  states 
that  the  form  of  the  equation  may  be  justified,  using  exact  many-body 
theory. 

To  determine  an  appropriate  form  for  f(q),  we  consider  the  calcula¬ 
tion  of  the  effective  interaction  between  quasi-  par  tides  in  an  inter¬ 
acting  electron  gas.  As  discussed  in  Ref.  8,  in  the  random-phase  ap¬ 
proximation  and  the  high-density  limit,  a  partial  matrix  element  summa¬ 
tion  of  "ring”  Feynman  diagrams  is  appropriate.  The  effective  inter¬ 
action  potential  is  calculated  to  be 


Veff  q 


4iTe 


2 


2 

q 


(66) 


which  in  real  space  is  a  screened  coulomb  interaction  in  the  Thomas- 

Q 

Fermi  approximation, 


2  -kr 

VeffODocf-e  3 


o  in 

Sham  follows  Hubbard’s  treatment  and  replaces  ^  by 


4ire 


(67) 


V  2  2  2  * 

eff  q  q2  +  kj  +  kg 


(68) 


where  kf  is  some  average  of  the  Fermi  vector.  Hubbard  also  proposes 

that  the  effect  of  exchange  on  the  screening  in  the  high-q  limit  should 

be  to  cancel  half  the  direct  coulomb  contribution  to  the  screening.  We 
2,7,10,11 


i 

4 


% 

■Vt 


will  write 
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eff  q 


1  /  4ire2  \ 

"  2  l  2  _  2  ’ 

+  9W 


2  2  2  2 
where  we  have  replaced  k^  +  kg  by  £k^.  This  is  in  accord  with  Wallace. 

5  is  an  adjustable  parameter  that  we  will  determine  later  to  ensure  that 

2 

I (q)  has  the  correct  q  +  0  limit. 

The  final  form  for  f (q)  to  be  used  in  Eq.  (65)  is,  therefore. 


f  (q)  *  J  “2 


q  +  fr* 


We  have  specified  the  terms  in  the  energy  wave-number  characteris¬ 
tic  [Eq.  (65)],  including  exchange  and  correlation  effects  within  our 
level  of  approximation.  Heine^  discusses  the  exchange  and  correlation 
effects  and  states  that  the  uncertainties  arising  from  the  approxima¬ 
tions  used  in  determining  Eq.  (70)  are  greater  than  the  uncertainties  in 
the  pseudo potential .  We  continue  by  choosing  an  appropriate  local 
bare-ion  potential  Wfi^  to  complete  our  determination  of  F(q)  and 

6.  A  Model  Pseudopotential.  We  approximate  the  local  bare-ion 


pseudopotential,  wfi(r),  by  a  simple  model  which  behaves  like  the  poten- 

2 

tials  in  the  real  metal,  w  (r)  is  composed  of  two  contributions, 

o 


w„(r)  «  w  (r)  +  w  (r)  . 

8  z  c 


w  (r)  is  due  to  the  coulomb  attraction  of  the  ion  and  has  the  form 
z 


/  \  -ze 

w  (r)  »  - 

z  r 


with  Fourier  transform  (i.e.,  the  w  matrix  element) 

zq 


k 


i 

i 

i 

i 


t 


n 


m  ^TfZe 


zq 
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(72) 


w^(r)  Is  a  localized  potential  confined  to  the  core  region.  It 
approximates  the  repulsion  due  to  the  core  electrons  and  tends  to  cancel 
the  coulomb  part  within  the  core.  An  appropriate  choice  for  wc(r)  Is  of 


the  form  of  the  s-state  core  functions. 


W  (?)  -  e-6r 


The  Fourier  transform  of  this  is 


1.2 


w 


,  Sail 

cq  % 


l<62  +  ,l>1] 


We  introduce  the  arbitrary  constants  g  and  p  such  that 


cq  0. 


Y) 


[(1  +  q2p2)'] 


(73) 


The  form  that  we  will  use  for  the  matrix  element  of  the  bare  ion  po¬ 
tential  in  Eq.  (65)  is,  therefore. 


„  -  _L  r-‘W2  +  — & — _i 
Bq  ao  L  q2  (l  +  q2p2)  J 


(74) 


7.  Total  Ion- Ion  Interaction.  For  a  molecular  dynamics  calcula¬ 
tion  we  are  interested  in  the  total  ion- ion  interaction  potential.  We 
have  just  calculated  the  effective  interaction  clue  to  the  presence 

of  the  nearly  free  electron  gas  in  a  simple  metal.  To  this  we  must  add 
the  coulomb  replus ion  between  the  lons,^ 


$ 


s 

* 


* 

I 

* 


£ 


and  the  exchange  repulsion  betveen  two  ion  cores 


i 


■j 


V 


-v 


(76) 


2 

which  is  a  Born-Mayer  repulsive  central  potential.  ot^  and  are 
empirically  determined  constants. 

Therefore,  we  may  summarize  the  results  of  this  development  by 
writing  down  the  ion- ion  interaction  potential,  which  will  be  used  in  our 
molecular  dynamics  calculations : 


<t>(  r) 


22  ^Br 

r  +  v 


(77) 


where 


w>  -  jC° F<q>  q2  d<* 


(78) 


F(q) 


e(q)  -  1  ■ 


2  *  Bq  1  +  [e(q)  -  1][1  -  f(q)]  ’ 


.  V 

i„  1 2  . 

87re2 

'  Bq  •  : 

2 

/,  2 

me  i 

ft  -  n 

2  2' 
2irkfft  rf 

(  2n 

g(q)  ~  * 


In 


1  t-a 


l  -  n 


+  l)  * 


2k, 


(79) 


(80) 


f  (q) 


2(q2  +  £k*) 


and 


1  f-Airze^  .  8  1 


(81) 


(82) 
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The  determination  of  the  constants  otB,  £,  6 ,  and  p  will  be  discussed 
in  Sec*  III. 

B.  Molecular  Dynamics 

We  have  calculated  an  effective  interaction  potential  for  the  ions 
in  a  simple  metal  [Eq.  (77)] .  We  use  this  as  the  potential  of  inter¬ 
action  between  pairs  of  classical  particles  and  will  solve  for  the  ion 
motion  in  the  metal  using  the  molecular  dynamics  technique.  In  this 
section  we  present  the  equations  that  are  solved  by  high-speed  com¬ 
puters  to  yield  ion  positions  and  velocities.  We  derive  the  conserva¬ 
tion  of  energy  and  momentum  in  a  molecular  dynamics  system,  describe  sys¬ 
tem  size  and  boundary  effects,  and  discuss  the  methods  of  calculating 
thermodynamic  properties  from  the  molecular  dynamics  results. 

1.  Central  Difference  Equations.  For  a  three-dimensional  array 
of  N  identical  particles,  we  solve  the  set  of  classical  Newtonian  equa¬ 
tions, 

mr^t)  «  fjU)  ,  (83) 

where 

i  -  particle  number  -  1,2,-**,N  , 

-  particle  position  (x^y^z^  , 

m  *  particle  mass,  and 

»  force  on  the  ith  particle. 

The  dots  denote  time  differentiation.  The  force  on  the  ith  particle  is 
determined  generally  by  the  positions  of  the  other  particles,  and  thus 
the  equations  are  coupled.  We  develop  the  difference  equations  used  in 
the  molecular  dynamics  program  to  solve  this  system  of  equations.  Con¬ 
sider,  for  simplicity,  only  the  x-direction.  The  y-  and  z-direction 
equations  are  identical. 
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In  the  x-direction 


"«i(t)  -  Fxl(t)  •  (84) 

Performing  a  Taylor  series  expansion  on  At,  the  new  position  at  t  +  At 
is 

xA(t  +  At)  -  x±(t)  +  At  X^t)  +j(At)2x1(t)  +  •••  ,  (85) 


and  similarly, 

x^t  -  At)  -  x^t)  -  At  x^t;  +  y(At)2x  (t)  +  •••  . 
Adding  these  equations  yields 

x^t  +  At)  -  -xi(t  -  At)  +  2xi(t)  +  (At)2x1(t)  +  •••  . 
Subtracting  gives 


x±(t) 


x^t  +  At)  -  xi(t  -  At) 
2At 


(86) 


(87) 


(88) 


These  equations  constitute  a  straightforward  central  difference 
scheme  that  is  appropriate  to  the  solution  of  the  i.th  particle  position 

3 

and  velocity  and  is  accurate  to  the  order  of  (At)  . 

We  work  easily  in  terms  of  displacements  defined  at  the  time  Inter¬ 
val  midpoints  as 

Axt(t  +  4p)  -  xt(t  +  At)  -  x^t)  .  (89) 

With  this  definition,  and  expanding  x.(t  +  At/2)  in  a  Taylor  series 
about  At/2,  the  following  difference  equations  are  calculated  as 
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2  Fxi(t) 


A*i(t  +  ^r>  ■  -  t5  +  (At)‘ 


xi(t  +  At)  -  xi(t)  +  Axi(t  +  4p)  ,  and 


Xi(t  +  T> 


At  +  (At/2)] 


At 


(90) 

(91) 

(92) 


where  we  have  used  x^(t)  ■ 

These  are  the  equations  used  in  the  molecular  dynamics  program  to 
solve  the  coupled  system  of  differential  equations.  The  force  of  the 
ith  particle  is  calculated  by  assuming  additive  interactions  between 
pairs  between  the  irh  particle  and  its  neighbors. 


P*l<‘> 


(93) 


The  prime  on  the  summation  indicates  that  the  j»i  term  is  omitted. 

is  the  x-component  of  the  force  on  the  ith  particle  due  to  the  ith 
particle  and 


4> 


F  -  -F 
xij  xji 


(94) 


by  Newton’s  third  law.  For  this  study  we  calculate  the  force  by  taking 
the  negative  gradient  of  the  ion-ion  interaction  potential  determined 
through  pseudopotential  theory  [Eq.  (77)]. 

The  computer  program  that  solves  the  difference  equations  is  easily 
understood .*  After  determining  the  initial  positions  and  velocities  of 
the  particles  (see  Sec.  IV),  all  the  particle  forces  are  calculated 

* 

The  molecular  dynamics  computer  program  used  in  this  study  was  developed 
by  B.  L.  Hollan  and  G.  K.  Straub  of  the  Los  Alamos  National  Laboratory. 
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baaed  on  the  current  positions.  With  Fxi(t)  thus  calculated, 

Ax^(t  +  At/2),  x^t  +  At),  and  x^t  +  At/2)  are  evaluated  readily  by 

Eqs.  90-92.  With  x^t  +  At)  known,  the  time  step  is  advanced  and 

Fxi(t  +  At)  is  calculated.  In  this  manner,  the  trajectories  of  all  N- 

3 

particles  are  calculated  exactly  at  any  time  [within  (At)  ].  The  par¬ 
ticle  positions  and  velocities  may  be  evaluated  to  obtain  a  time  history 
of  the  system  being  studied. 

2.  Conservation  of  Momentum  and  Energy.  The  difference  equations 
just  described  can  be  shown  to  explicitly  conserve  momentum  and  energy.* 
In  a  system  of  N  identical  particles  interacting  by  way  of  a  conserva¬ 
tive,  additive  force  between  pairs,  where 

-34>(x  ) 

Fxi(t)  •  “i(t)  *  Tx7  -  (95) 

and  4>  is  the  additive  potential  between  pairs,  we  can  demonstrate  con¬ 
servation  of  momentum,  P,  by  showing  that 


Ht  *  f) 


or 


SI ■ii(t  ■r  T>  "H  “i(t  "  *£> 

i-l  i-1 


(97) 


Using  the  difference  equations  [Eqs.  (90)  —  (9 2) ]  we  evaluate 


N 


21  V* 

i-l 


*f> 


+  t> 


-  f  >  + 


Based  on  a  calculation  by  Brad  Lee  Hollan  of  the  Los  Alamos  National 
Laboratory . 
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Also, 


N  N  N  N  N  . 

£  r*i  -££  F„ij  -  EE  <Fxu  +  v* 

i-1  i-1  j-1  i-1  j>l 


But,  by  Newton’s  third  law 


F  *  F 
xij  xji  ’ 


so  that  this  sum  is  zero  and  we  have  shown  in  Eq.  (97)  that 

£  mxi(t  +  ^r)  tnxi(t  -  )  , 

i-1  i-1 

which  proves  that  momentum  is  conserved  to  within  the  accuracy  of  the 
difference  equations. 

Similarly,  we  calculate  the  total  energy  to  demonstrate  conserva¬ 
tion  of  energy. 

E  -  total  energy  -  potential  (PE)  +  kinetic  (KE)  energy. 

We  calculate  the  change  in  potential  energy, 


APE  -  PE(t  +  4p)  -  PE(t  -  jf- )  , 


and  compare  it  with  the  change  in  kinetic  energy. 


ARE  -  KE(t  +  yO  -  KE(t  -  y)  . 


The  change  in  potential  energy  is 


j  rtf’ll n'l' f  i~'"  A  '  “A  •'  4* Vlkxtimt&ltMm*  '■■ 


V 


We  evaluate  these  terms  by  noting  that 


+  ^)2  -  xt( t  -  -  U1(e  +  +  x1(t  -  ^)J 


I*t(t  +  yO  -  x1(t  -  ^)] 


and  by  expanding  using  the  Taylor  series 


2  t 

X±(t  +  Y')  -  xt(t)  +  Y  **(*)  +  j(^)  ’^(t)  +  O(At)' 


x^t  -  Y>  -  x^t)  -  X^t)  +  j(^)2  x1(t)  +  O(At)3 


By  adding  and  subtracting  the  second  of  these  equations  from  the  first* 
we  get  the  needed  factors  in  Eq.  (102).  By  using 


F±(t)  -  23ixi 


and  some  algebra,  we  arrive  at 


AKE  «>  [Fi(t)xi(t)  +  0(At)  j  , 


and  we  see,  from  Eqs .  101  and  103,  that 


AE  -  AKE  +  APE  -  0  +  0(At)J  . 


Thus  we  have  shown  explicitly  that  both  momentum  and  energy  are  con¬ 
served  by  the  difference  equations. 

3.  Boundary  Conditions  and  System  Size.  We  now  set  up  the  system 


of  N  interacting  particles.  We  perform  our  molecular  dynamics  calcula¬ 
tion  on  this  system,  and  from  the  resulting  information,  we  infer  the 
macroscopic  properties  of  sodium.  The  system  is  obviously  limited  to 
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having  N  much  less  than  a  macroscopic  (10  )  number*  We  must  minimize 

the  effects  of  finite  system  size  by  choosing  appropriate  boundary  con¬ 
ditions,  and  we  must  recognize  the  N-dependence  of  our  results. 

To  minimize  the  effects  of  finite  system  size,  we  choose  periodic 
boundary  conditions  whereby  the  N-particle  system  is  repeated  periodic¬ 
ally  throughout  space  and  thus  may  be  considered  an  infinite  system  with 

the  imposed  periodicity.  Figure  1  is  a  schematic  of  an  N*3  system  in 

12 

two  dimensions  and  its  nearest  repeated  systems.  L  is  the  length  of 
the  square  system  box  and  RMAX  is  the  range  of  the  potential.  We 
use  the  minimum  image  convention,  which  means  that  a  particle  only 
interacts  with  the  image  of  a  neighbor  that  is  nearest  to  it.  That  is,  in 
Fig.  1,  particle  two  will  only  interact  with  the  image  of  particle  three 
in  the  box  to  its  right  and  no  other  "particle  threes.1*  This  constrains 
us  to  keep  the  system  length,  L,  in  a  given  direction  greater  than  twice 
the  potential  range  so  that  not  only  will  a  particle  never  interact  with 
its  own  image  but  it  will  not  interact  with  more  than  one  of  its  neigh¬ 
bor’s  images  *  The  method  for  determining  our  particular  system  dimen¬ 
sions  and  the  potential  range  is  discussed  in  Sec.  IV. 

The  N  dependence  of  calculated  system  properties  has  only  been 
determined  for  some  simple  cases.  (See  discussions  in  Refs.  12,  13, 

and  14.)  Based  on  such  studies,  we  expect  a  1/N  dependence.  However, 

12 

a  generalized  solution  method  does  not  exist.  Therefore,  we  must  in¬ 
clude  sufficient  particles  in  the  system,  and  we  check  to  see  that  the 
calculated  results  are  not  N  dependent.  If  there  is  an  N  dependence 
we  may  be  able  to  evaluate  how  to  extrapolate  our  results  to  the  macro¬ 
scopic  N  limit.  In  our  equilibrium  property  studies,  we  are  not  overly 
concerned  because  the  N  dependence  is  generally  small  for  systems  with 
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Fig.  1. 

A  two-dimensional  system  of  N*3  particles  in 
an  L  x  L  box  with  periodic  boundary  conditions. 
RMAX  is  the  range  of  the  potential.  Note  that 
as  particle(T)  moves  out  of  the  system  it  re¬ 
appears  by  way  of  its  periodic  image. 


7+ 


7* 
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n  >  100  and  we  will  be  dealing  with  N  >  600,  owing  to  the  long  range  of 
the  pseudopotential.  If  we  study  nonequilibrium  properties  or  effects 
relating  to  large  numbers  of  particles,  such  as  coexistence  of  phases, 
then  the  N  dependence  becomes  more  critical.  We  discuss  N  dependence 
effects  as  they  affect  our  calculations  in  Sec.  IV. 

4,  Thermodynamic  Properties 

When  it  is  appropriate  to  calculate  macroscopic  equilibrium  prop¬ 
erties  from  an  N-particle  molecular  dynamics  system,  we  proceed  by  tak- 

12 

ing  averages  over  time.  In  a  molecular  dynamics  calculation,  the 
number  N  of  particles,  volume  ft  of  the  system,  and  total  energy  E  remain 
constant.  With  these  constraints,  the  system  forms  a  microcmonical  en¬ 
semble.  However,  there  is  the  additional  constraint  of  constant  linear 
momentum  (M) ,  so  we  would  describe  the  system  as  having  N,  V,  E,  and  M 
constant. ^  This  must  be  taken  into  account  when  comparing  molecular 
dynamics  time  averages  with  ensemble  averages  obtained  using  statistical 
mechanics  or  other  calculational  techniques  such  as  Monte  Carlo  (see 
Refs.  12,  13,  and  15  for  discussion).  We  are  not  concerned  immediately 
with  such  comparisons.  The  essential  assumption  regarding  the  appropri¬ 
ateness  of  the  averages  is  that  the  system  is  ergodic.  That  is,  all 

13 

states  of  the  system  in  phase  space  are  mutually  accessible  with  equal 
probability,  so  that  all  time  histories  of  the  same  system  have  equiva¬ 
lent  statistical  averages. ^  This  assumption  has  never  been  proved, 

12 

but  experience  shows  it  to  be  valid.  We  assume  that  our  time  averages 
yield  the  appropriate  statistical  averages  of  the  system’s  equilibrium 
properties. 

We  will  now  determine  various  properties  of  our  system  at  a  partic- 


* 


ular  time,  with  the  understanding  that  we  will  sample  these  properties 
many  times  during  a  molecular  dynamics  calculation  and  average  them. 


We  calculate  the  average  kinetic  energy  per  particle  of  the  system 


of  N  particles  as 


‘s  E  I  m(*2i  +  +  Z1}  • 


The  average  potential  energy  per  particle  is 

i-i  j-i 


where 


ri  "  (xi»y1»zi) 


and  ((>(r)  is  the  potential  that  describes  the  interaction  between  two 
particles.  The  factor  of  1/2  is  included  to  correct  for  double- summing 
of  the  potential  energies. 

Our  system  at  equilibrium,  obeying  the  laws  of  classical  dynamics, 
will  satisfy  the  Maxwell-Boltzmann  law.  The  law  states  that  the  average 
of  a  property  (P,  for  example)  can  be  obtained  by 


-  f  MS  l  dV"JP3N  , 

8  'n-j*-,,kT‘v<p 3» 


where  E  is  the  system  energy,  kT  is  Boltzmann's  constant  times  the  tem¬ 
perature,  and  the  integration  is  over  all  points  in  phase  space  with 
positional  coordinates  qi  and  momentum  coordinates  p^.^ 

Using  this  law,  we  may  calculate  the  distribution  function  for 


each  velocity  component  as 


*  A.-  “*v  ■  1  v.  :  i.  k  > 


IN.  .  .  .2, 
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and  similarly  for  y  and  2.  This  is  the  Maxwell  or  normal  distribu¬ 
tion.*^  Then  we  calculate  the  average  kinetic  energy  per  particle  as 


KE  -  |(x2  +  y2  +  z2)  -  |  kT  . 


We  are  assuming  that  our  system  has  no  net  motion,  i.e.  x  *  y  «  2  *  0. 
If  there  were  net  motion,  we  would  calculate  only  the  fluctuations 

2  t  2 

about  this  average  motion  and  replace  x  by  (x  -  x)  .  This  connects 
the  kinetic  energy  per  particle  calculated  using  molecular  dynamics  and 
the  thermodynamic  temperature  of  the  system. 


1  V  1  ,.2  ^  .2  ^  .2.  3 

n  2^  2  m(xi  +  y±  +  zi)  *7 


kT  . 


We  also  accept  this  expression  as  the  definition  of  kinetic  temperature 

18 

to  be  used  even  when  the  system  is  not  in  equilibrium. 

Note  that  the  N-particle  molecular  dynamic  system  may  not  be  in 
equilibrium.  Also,  the  velocity  distribution  may  net  be  normal.  We  may 
investigate  this,  as  far  as  the  temperature  and  kinetic  energy  are  con¬ 
cerned,  by  seeing  if  the  kinetic  temperature  is  isotropic.  That  is, 
we  may  for  convenience  define 


1  ,T  1  A  1  .2 

2  kTx  "  N  L  2  “1  * 
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and,  similarly,  T  and  T  .  We  require  that  T  -  T  -  T  for  the  system 
y  2  x  y  z 

to  be  in  equilibrium.  To  investigate  if  the  velocity  distribution  is 
normal,  we  may  calculate  the  kurtosis  (which  measures  the  "degree  of 

19 

peaking")  of  the  distribution  in  the  molecular  dynamics  calculation. 

We  calculate  the  kinetic  temperatures  and  the  kurtosis  for  example 
calculations  and  discuss  their  implications  in  Sec.  IV. 

The  volume  of  a  molecular  dynamics  calculation  as  we  apply  it  is 
a  constant  and  an  input  parameter.  A  given  calculation  will  yield  the 
potential  and  kinetic  energies  and  the  calculated  temperature.  For 
each  calculation  we  have  an  equation-of -state  point  for  the  molecular 
dynamics  system  consisting  of  a  total  internal  energy  (kinetic  plus 
potential)  at  a  volume  and  temperature. 


Hr.  A.w  V  >  -V  v-4  •„  ~  > i  ;ht  tyMMA  ? 


V, 
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III.  THE  ION-ION  INTERACTION  POTENTIAL  FOR  SODIUM 

In  this  section  we  discuss  how  the  necessary  parameters  are  deter¬ 
mined  for  the  calculation  of  the  ion- ion  interaction  potential  in  Eq. 
(77).  We  use  these  parameters  to  calculate  the  functions  in  Eqs.  (78)  to 
(82)  which  specify  the  potential.  Finally,  we  calculate  the  interaction 
potential  and  discuss  how  it  is  tabulated  for  use  by  the  molecular  dy¬ 
namics  program. 

A.  Determination  of  the  Parameters 

1.  Determination  of  £.  The  constant  £  in  the  expression  for  the 
function  f(q)  [Eq.  (81)],  which  corrects  the  dielectric  function  to  in¬ 
clude  exchange  and  correlation  effects,  may  be  determined  analytically. 
In  the  long  wavelength  (q  +  0,  or  macroscopic)  limit,  the  static  di¬ 
electric  constant  is  related  to  the  electron  gas  compressibility,  K, 


£ 1  (q) 


,  xt2  2 

1 

q 


q  0  limit 


(111) 


where  the  prime  distinguishes  this  dielectric  constant  from  the  Hartree 
dielectric  constant  of  the  noninteracting  electron  gas,  £(q),  defined 
by  Eq.  (80). 

Note  that 


Zjf.  -  !L- 

(e'(q)  -  1)  k 


(112) 


where  the  superscript  zero  indicates  the  values  appropriate  to  a  non¬ 
interacting  electron  gas.  With  some  algebra  and  the  formulas  relating 
f (q)  and  the  dielectric  constants  [for  example,  Eq.  (63)],  this  relation 
becomes 


I 


~  -  1  -  f(q)[£'(q)  -  1]°  -  1  -  f(q)[£(q)  -  D  , 


which  is  valid  in  the  long  wavelength  limit. 

From  the  expression  for  e(q),  Eq.  (80),  we  find** 


4me2k~ 

lim  [ e(q)  -  1]  -  — YY 
q-*0  Trfi  q 


3irh  2k, 


so  that 


-  lim.  [1  -  f(q)(£(q)  -  1)] 
q-0 


2 

4me  kr  2 

f  me 


2 (q2  +  £k2)/  WV  37r!i2kf 


,  2me2 

K  £kf  7T?> 2 


This  is  the  relation  that  determines  £. 

We  calculate  <  and  from  the  relation 


i 

<-*0  727  •  (117) 

a 

where  EQ  is  the  electron  gas  ground-state  energy.  We  write  the  energies 

in  terms  of  r  ,  which  is  defined  as  the  radius  of  the  sphere  whose 
s 

volume  is  the  average  conduction  electron  volume,  so  that 


4  3  0 

■=■  Ttr  *  — 
3  s  z 
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To  calculate  the  compressibilities,  we  use  the  results  from  Sec.  V 
for  the  ground-state  energies  of  the  noninteracting  and  interacting 
electron  gas.  Eg  and  EQ,  respectively. 


E°  -  2‘21 
E0  r  2 
s 


Ry 


-  -  °-;9-16  -  0.115  +  0.031  la  r  \R 
.  r2  rs  sl  y 


With  the  help  of  the  relations 


we  may  calculate  <°/<  and  use  Eq. (116)  to  arrive  at 

5  -  0,916/(0,458  +  0.012  r  )  .  (119) 

3 

With  the  value  -  3.939  a^  corresponding  to  an  atomic  volume  of 

3 

256  a^  for  sodium,  we  have 

5  -  1.81  .  (120) 


We  notice  here  that  $  is  volume-dependent  and  is  treated  as  such  in  our 
calculations . 

2.  Determination  of  3,  and  p.  The  remaining  parameters 

Ygi  S,  and  p  have  been  determined  by  fitting  to  available  data. 
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We  have  chosen  the  value  for  in  the  Born-Mayer  repulsive  poten¬ 
tial,  Eq.  (76),  to  have  the  value  determined  by  simultaneously  fitting 

calculated  lattice  energies  to  equation-of -state  data  for  the  family  of 
20 

alkali  halide  salts.  The  value  is 


—  -  0.339  x  10-8  cm 


so  that 


y 


B 


1.56  a 


-1 


(121) 
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Wallace  determined  the  values  of  otg,  B>  and  p  by  fitting  the  cal¬ 
culated  expressions  for  the  total  adiabatic  potential  and  its  volume 
derivative  to  measurements  of  the  equation-of -state  properties  of  sodium. 
The  data  included  were  the  binding  energy,  the  ionization  energy,  and 
the  bulk  modulus  at  zero  temperature  and  pressure.  The  requirement  that 
the  pressure  be  zero  at  zero  temperature  was  also  used.  He  found  that 
these  data  could  be  fit  with  some  arbitrariness  remaining.  This  arbi¬ 
trariness  was  removed  by  requiring  that  the  calculated  average  of  the 

2 

phonon  frequencies  squared  (<o>  >  as  calculated  by  lattice  dynamics) 
also  be  fit.  The  resulting  parameters  are 


ctB  *  10.5  Ry  , 

3 

8  *  37.  Ry  Sq  ,  and 
p  -  0.50  a^  . 

We  discuss  in  detail  the  total  system  energy  calculation  and  equa- 


tion-of-state  value  determination  in  Sec.  V. 
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With  these  parameters  determined,  we  have  completely  specified  the 
effective  ion-ion  interaction  potential  of  Eq.(77).  We  now  proceed  to 
calculate  this  potential  and  apply  it  to  the  molecular  dynamics  calcula¬ 
tions. 

B.  Calculation  of  the  Potential 

We  now  calculate  the  ion-ion  interaction  potential  for  sodium,  using 
Eqs.  (77)-(82)  and  the  parameter  values  determined  in  the  previous  sec¬ 
tion.  For  these  calculations,  which  illustrate  the  factors  involved  in 
the  ion- ion  interaction,  we  use  the  observed  zero  temperature  and  pres¬ 
sure  atomic  volume  for  sodium  of 

Q0  »  256  Sq  , 

which  yields  a  Fermi  wave  vector  of 

/  2  \l/3 

kf  ■  pif)  -  °-4872  • 


1.  e(q).  The  static  Hartree  dielectric  function  is  calculated  us¬ 

ing  Eq.  (80)  and  the  result  is  the  solid  line  shown  in  Fig.  2.  It  acts 
as  the  dielectric  function  for  an  interacting  electron  gas  without  tak¬ 
ing  exchange  and  correlation  into  effect.  It  has  large  and  small  q 
2 

limits  given  by 


and 


(123) 


(124) 


Fig.  2. 

The  solid  line  is  the  static  Hartree  dielectric  function,  e(q), 
as  given  by  Eq.  (80).  The  dashed  line  is  the  modified  dielectric 
function  as  given  in  the  discussion  following  Eq.  (125). 


e(q)  has  a  logarithmic  infinity  in  its  second  derivative  at  q  ■  2 
which  will  have  an  important  effect  on  the  interaction  potential,  as 
will  be  discussed  later. 

2.  f (q) .  The  interpolation  formula  f(q),  which  corrects  the 
electron  screening  because  of  exchange  and  correlation  effect,  is  given 
by  Eq.  (81)  and  is  shown  in  Fig.  3.  It  varies  smoothly  from  the  expected 
f(0)  *  0  long-range  limit  to  the  f(°°)  *1/2  limit,  where  the  screening 
effect  is  reduced  by  a  factor  of  2  because  only  electrons  of  anti¬ 
parallel  spin  should  interact  in  this  limit  with  exchange  taken  into 
account. ^  With  f(q)  included,  the  relation  between  the  total  potential 
W  and  the  bare- ion  potential  W_  Fourier  coefficients  is 

o 

wn  *  +  <c<q>  -  1X1  -  f(q)]  >  (125) 

q  aq 

so  that  the  modified  dielectric  function  is  1  +  [e(q)  -  1]  [1  -  f(q)]. 
This  function  is  shown  as  the  dashed  line  in  Fig.  2. 

3.  wD  and  w  .  The  local  bare-ion  pseudopotential  matrix  element 

_ _ 1 

w  ,  discussed  in  Sec.  I  and  given  by  Eq.  (82)  is  the  solid  line  in  Fig. 
bq 

4.  The  screened  matrix  element  given  by  Eq.  (125)  is  the  dashed  line  in 
the  same  figure. 

Note  here  that  the  screening  effect  is  the  cutting  off  of  the  long- 
range  (small  q)  part  of  the  potential,  as  expected.  We  also  note  that 
in  the  small  q-limit 


w  -  =i 

q  T\2m  *f/  3 


ef  s  -0.158  Rv  , 


where  we  have  used  the  liiji  f(q)  ■  0  and  the  small  q-limit  for  the 
Hartree  dielectric  function  given  by  Eq.  124. 


Fig.  3. 

iterpolation  formula  for  the  approximate  correction  to  the  %; 

on  screening  due  to  exchange  and  correlation  effects9  f(q), 

en  by  Eq.  (81),  jfc 

* 


* 


or  w  (Ry) 


f 


I 

I 
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Fig.  4. 

The  solid  line  is  the  local  bare-ion  pseudopotential  matrix  ele¬ 
ment,  wgq,  given  by  Eq.  (82).  The  dashed  line  is  the  screened 
matrix  element  given  by  Eq.  (125). 
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4 .  F(q) .  The  energy  wave-number  characteristic  given  by  Eq.  (79) 
is  calculated  easily  using  the  results  above,  and  it  is  plotted  as  a 
function  of  q  in  Fig.  5.  F(q)  has  more  detail  than  this  scale 
graph  shows,  such  as  a  second  negative  "hump"  at  about  q  *  Sa^  with  a 
magnitude  of  1.0  x  10  7  Ry. 

In  the  large  q-limit,  F(q)  goes  to  zero  as 


lim  F(q) 

q-*oo 


2  3 
64  z  k^ 


and  in  the  small  q-limit,  F(q)  goes  to  negative  infinity  as 


0  2  2 
-2ttz  e 


5.  V-%TTV(r).  We  can  evaluate  the  effective  ion-ion  interaction 
_ IND 

due  to  the  presence  of  the  electrons  in  the  ion-electron  system.  We 
must  evaluate  the  integral  in  Eq.  (78),  which  is 


V„m(r) 


^0  ^  sin  qr  2  , 

^JQ  • 


It  is  instructive  to  look  at  plots  of  the  argument  of  this  integral; 


AEG  -  '-f  F(„)  ,2 

TT  “  ** 


which  are  plotted  for  different  r  values  in  Fig.  6.  This  figure  shows 
that  the  absolute  value  of  the  integral  will  decrease  with  increasing 
r  because  of  the  modulation  of  the  sin  (qr)  term. 
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Fig.  6. 

Arguments  of  the  effective  ion-ion  interaction  integral  as  given 
by  Eq.  (130)  plotted  for  different  r  values. 


V 

1 


The  integral  of  this  argument  is  evaluated  numerically  with  the 
help  of  a  Simpson's  rule  integration  subroutine  and  the  following  re¬ 
lationship. 


Jr00  r  b  /• 00 

f  ARG  dq  -  /  ARG  dq  +  /  (lim  ARG)  dq  , 


where  the  second  integral  is  an  analytical  integral  of  the  large  q-limit 
of  the  argument. 

Using  the  large  q-limit  for  F(q)  given  by  Eq.  (127),  we  must  evaluate 


(lim  ARG)  dq 
q-H» 


,,  3 

-64  z  k^r 


sin  x 


'br  x 


where  we  have  changed  variables  in  the  integral  from  qr  to  x.  This 
integral  may  be  evaluated  using  standard  integration  formulas  and  the 
rational  approximations  to  the  sine  integral  of  the  form 


sin  x 


given  by  Ref.  22. 

In  practice,  the  first  integral  in  Eq.  (131)  is  broken  into  NINT 
intervals,  starting  at  A,  with  each  B  wide  to  ensure  that  the  Simpson's 
rule  integration  is  able  to  converge  efficiently.  The  integrator  uses 
a  convergence  criterion  parameter,  EPS,  to  determine  the  accuracy  of 
each  integration.  The  actual  integration  procedure  is  illustrated  by 


c  oo  i*  A+B  ^»A+2B 

J  v  +J 

•'n  yA  J  a 


0  'A 


rA+2B  /-A+(NINT)B  r  ® 

+  •*•/  +  j  (large  q  limit)  .  (133) 

A+B  •'A+CNINT-DB  /a+(NINT)B 
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Breaking  up  the  interval  allows  investigation  and  correction  for  in¬ 
accuracies  that  may  arise  when  Simpson’s  rule  is  used  to  evaluate 
integrals  that  have  positive  and  negative  areas  canceling.  We  performed 
a  convergence  study  of  the  V^p(r)  integral  and  determined  that  the 
values 


NINT  -  20  , 

B  -  1  , 

A  ■  1  x  10  ^  ,  and 


EPS  *  1  x  10 


yielded  results  accurate  to  about  10  »  relative  to  the  integral  values 

we  obtained  for  the  well-converged  solution.  This  value  set  was  chosen 
to  evaluate  the  integrals  for  use  in  the  rest  of  the  study.  The  func¬ 
tion  V_^(r)  calculated  in  this  manner  is  shown  in  Fig.  7  and  as  the 
INI) 

solid  line  in  Fig.  8, 

We  now  investigate  the  Vj^(r)  behavior  in  the  large  r-limit.* 

To  do  this  we  expand  Eq.  (79)  forF(q)  by  letting 


_ i(q)  -  i _ - _ * _ 

1  +  [G (q)  -  1][1  -  f (q) ]  1  -  f(q)x  ’ 

where 


e(q) 


We  note  from  the  values  of  f(q)  and  E(q)  that  0  <  f(q)  <1/2  and  0  < 
x  <  1,  so  that  [f(q)x]  <  1.  We,  therefore,  may  expand  Eq.  (135)  as 


& 

Based  on  a  calculation  by  Galen  K.  Straub,  Los  Alamos  National 
Laboratory. 


as  given  by  Eq.  (78). 


Fig.  7. 

The  effective  ion- ion  interaction, 


{ 

9 

1 


Fig.  8. 

Plots  of  the  effective  ion- ion  interaction,  (solid  line); 

the  coulomb  term,  2^e2/r  (dashed  line);  the  Born-Mayer  term,  age~^®r 
(dotted  line);  and  their  sum,  the  total  ion-ion  interaction,  <p(r) 
(chain-dashed  line)  as  given  by  Eq.  (77). 
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The  first  term  corresponds  to  the  coulomb  contribution  of  the  bare 
ion  potential  and  cancels  the  direct  coulomb  ion-ion  interaction  at 
large  r.  The  second  term  represents  the  Friedel  oscillations  typical 
of  screening  caused  by  electrons  in  a  system  characterized  by  a  sharp 
cutoff  of  momentum  at  the  Fermi  surface. ^  This  term  gives  rise  to 
the  long-range  oscillatory  behavior  of  the  ion- ion  interaction. 

6.  $(r).  We  calculate  the  total  effective  ion- ion  interaction 

given  by  Eq.  (77)  as 

2  2  -Y  r 

<Kr)  -  +  «Be  B  +  vIND<r)  •  (140) 

The  coulomb  and  Born-Mayer  repulsive  terms  are  plotted  as  the  dashed 
and  dotted  lines,  respectively,  in  Fig.  8.  aac*  are  a^so 

plotted  in  this  figure.  The  nature  of  (r )  is  obscure  in  this  figure, 
other  than  the  cancellation  of  the  coulomb  repulsion  by  the  leading  1/r 
term  in  the  large  r  expansion  of  vrNp(r)  given  by  Eq.  (139). 

The  nature  of  <Kr)  is  apparent  from  the  plot  in  Fig.  9.  We  notice 
the  dominant  repulsion  for  r  <  ^,5  aQ,  the  minimum  in  the  potential  at 
r  -  8ag,  and  the  long-range  oscillatory  behavior  as  discussed  above. 

For  convenience,  we  include  here  a  plot  of  the  force  between  pairs, 

F(r)  -  (141) 

in  Fig.  10.  The  force  is  the  value  used  directly  in  the  molecular  dy¬ 
namics  calculations  to  determine  ion  motion.  Values  for  both  the  po¬ 
tential  and  force  as  a  function  of  interatomic  distance,  r,  are  tabu¬ 
lated  in  Appendix  A. 

These,  then,  are  the  effective  ion-ion  interaction  potential  and 
force  that  will  be  used  in  our  molecular  dynamics  calculations.  We 
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notice  that  they  are  dependent  on  the  atomic  volume  through  the 
Fermi  wavenumber  k^,  so  that  new  functions  are  required  at  each  sodium 
density  calculated.  Because  the  calculation  of  the  potential  and  force 
functions  is  lengthy,  we  will  tabulate  the  results  and  interpolate  from 
the  tables  during  a  molecular  dynamics  calculation.  This  procedure  is 
discussed  in  the  next  section. 

C.  Tabulation  of  the  Potential  and  Force 

Because  each  molecular  dynamics  calculation  with  the  number  of 
particles  that  we  are  considering  is  inherently  costly,  we  want  to  make 
the  potential  and  force  information  available  in  a  manner  which  mini¬ 
mizes  the  number  of  operations  required  by  the  molecular  dynamics  pro¬ 
gram.  To  do  this  we  tabulate  the  potential  and  force  values  and  then 
look  up  the  values  as  needed. 

There  is  a  tradeoff  between  different  interpolation  schemes.  Some 
schemes  are  accurate  for  a  relatively  small  number  of  points  in  the 
table,  yet  require  many  operations  to  obtain  a  value.  The  simplest  of 
schemes,  linear  interpolation,  requires  a  high  point  density  but  few 
calculations.  Because,  in  a  molecular  dynamics  program,  we  must  perform 
many  table  lookups,  and  storage  on  the  computer  we  are  using  (CDC  7600) 
is  not  a  serious  constraint,  we  have  decided  to  use  a  linear  interpola¬ 
tion  scheme  with  enough  points  in  the  tables  to  ensure  a  negligible 
loss  of  accuracy  owing  to  the  table  lookup. 

The  potential  and  force  functions  are  dependent  both  on  interatomic 
distance,  R,  and  atomic  volume,  V.  (For  this  discussion,  we  are  using 
notation  convenient  to  the  notation  used  in  the  computer  program.) 
Therefore,  we  must  set  up  our  tables  and  interpolate  in  both  the  R  and 


V  dimensions. 
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To  illustrate  table  setup  and  the  interpolation  scheme,  we  will 
outline  the  procedure  for  determining  the  potential,  POT,  and  force, 

F,  once  V  and  R  have  been  specified. 

The  volume  table  is  set  up  with  a  minimum  volume  VMIN,  a  maximum 
volume  VMAX,  and  a  certain  number  of  table  values  NV  (see  Fig.  11). 

The  constant  interval  between  values  is  given  by 

DELV  »  VM^  —  yMIN  .  (142) 

NV  -  1 


The  radius  table  is  set  up  in  a  similar  manner  with  RMIN,  RMAX,  NR,  and 
DELR  values  specified. 

V(IV)  is  the  volume  at  the  IVth  position  in  the  volume  table, 

R(IR)  is  the  radius  at  the  IRth  position  in  the  volume  table,  and 
POT(IR,IV)  is  the  potential  at  IR  and  IV,  which  are  calculated  from 
pseudopotential  theory  as  discussed  in  Sec.  III.B  [and  similarly  for 
F  ( IR ,  IV )  ]  . 

Now,  given  an  arbitrary  R  and  V  within  the  table  limits,  we  per¬ 
form  an  integer  divide  (i.e.,  truncate  the  division  to  an  integer)  to 
get 


IV 


V  -  VMIN  . 
DELV 


(143) 


so  that  we  know  that  V  is  positioned  between  V(TV)  and  V(IV  +  1)  in  the 
volume  table. 

At  each  of  these  volumes  we  find  IR  and  IR  +  1  and  interpolate  to 

find 

P0T1  -  POT  at  IV  and  R  , 


P0T2  •  POT  at  IV  +  1  and  R 


VMIN 


Vf^ 


VOLUME 


Fig.  11. 

A  schematic  showing  the  arrangement  of  the  volume-dependent  table. 
NV  is  the  number  of  points  in  the  table. 
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We  may  then  interpolate  on  the  volume  to  arrive  at  the  value  for  the 
potential  at  V  and  R  given  by 

V  *  v(IV3 

POT  -  P0T1  +  — ~  - 1  (P0T2  -  P0T1)  ,  (144) 

and  similarly  for  the  force 

F  -  F1  +  V"dexvIV)  (f2  ~  F1)  •  (145) 

We  determined  the  appropriate  density  of  points  for  these  tables. 
Calculations  of  the  potential  and  force  at  different  volumes  and  con¬ 
stant  radius  showed  an  almost  linear  relationship,  and  five  tables  were 
adequate  to  cover  a  range  from  10%  compression  to  10%  expansion  from 

normal  density.  For  the  R  tables,  a  value  of  NR  *  2000  was  chosen,  which 

-4 

yields  interpolated  values  that  are  precise  to  about  10  ,  relative  to 

calculated  values.  The  RMIN  and  RMAX  values  include  the  expected  inter¬ 
atomic  distances. 

The  following  values  were  chosen  for  the  table  setup. 


RMIN  -  4.0  aQ 
RMAX  -  30.0  aQ 
NR  «  2000 

3  (146) 

VMIN  *  230.01  a^  (10%  compression) 

VMAX  *  281.12  a^  (10%  expansion) 

NV  -  5 


Curves  obtained  from  the  tables  for  the  10%  compression  and  ex¬ 
pansion  volumes  are  shown  in  Fig.  12  to  illustrate  the  density  depend¬ 
ence  of  the  potential. 


001  0'S  O'O  0'4-  0*01-  0S1-  O'K- 

(XM)  (4)4, 


10%  Compression 
10%  Expansion 


4*0  9.0  8.0  10. 0  12.0  14.0  16.0  10.0  20.0  22.0  24.0  26. 0  20.0  30.0 

r(aQ) 

Fig.  12. 

The  total  effective  ion-ion  interaction  potential  for  the  10%  com¬ 
pression  (ft()  *  230  a^)  and  10%  expansion  (Qq  *  280  aj?)  conditions. 
The  oscillating  curves  at  large  r  have  been  multiplied  by  100* 
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IV.  SETUP  OF  THE  MOLECULAR  DYNAMICS  CALCULATIONS 

We  have  described  the  molecular  dynamics  calculations  and  the 
pseudopotential  method  by  which  we  calculate  the  effective  pair  inter¬ 
action  between  ions.  A  given  molecular  dynamics  calculation  proceeds 
as  described  in  Sec.  II. B,  with  the  forces  needed  in  Eq. (9C)  read  as  re¬ 
quired  from  the  tables  set  up  as  described  in  Sec.  III.C.  In  this  sec¬ 
tion  we  discuss  the  setup  of  our  particular  calculations  on  sodium, 
and  describe  the  units,  initial  conditions,  and  crystal  configurations 
for  the  hexagonal  close-packed  (hep)  and  body-centered  (bcc)  phases  of 
sod  iuxn . 

A.  Units 

We  specify  the  energy  (E 0) ,  distance  (X0) ,  and  mass  (M0)  units  for 
these  calculations  as 

E0  *  1  Rydberg  -  Ry  -  13.60559  eV  -  2.17971  x  10_il erg  -  e2/ 2aQ  , 

X0  -  1  Bohr  radius  -  aQ  -  0.529167  x  10~8  cm  »  0.529167  A  -  h2/me2  ,  (147) 

-23 

M0  -  1  molecular  weight  of  sodium  *  22.9898  g/mole  «  3.81731x10  g 

Note  that  in  these  units,  the  mass  of  a  sodium  ion  is  unity.  From  these 

specified  units  we  derive  the  time  and  velocity  units  as 

/M0Y5  -15 

t0  -  X0f— )  •  7.00281  x  10  s  , 

and  (148) 

v0  -  M  .  0.755650  cm/us  . 
t0 

3 

Pressure  is  given  in  Ry/a^.  We  will  use  kelvins  (K)  as  our  temperature 
unit  with  Boltzmann’s  constant  given  by 

k  -  6.33359  x  10-6  Ry/K  . 


■ 


;-V 
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B.  Initial  Conditions 

To  solve  the  difference  equations  described  in  Sec.  II. B,  we  must 
supply  initial  positions  and  velocities  for  the  particles  in  the  system, 
such  that  both  the  expected  crystal  configuration  and  approximate  tem¬ 
perature  are  predetermined.  In  the  next  section  we  discuss  the  deter¬ 
mination  of  the  initial  particle  velocities.  In  the  following  two 
sections  we  describe  the  initial  positions  for  hep  and  bcc  crystals  of 
sodium. 

1.  Particle  Velocities.  We  determine  the  initial  velocities  to 
satisfy  the  requirements  that  the  center  of  mass  velocity  is  zero  and 
that  the  average  of  the  velocities  squared  will  give  twice  the  tempera¬ 
ture  desired.  We  say  twice  here  because,  with  the  particles  placed  at 
their  lattice  positions,  all  Initial  energy  will  be  kinetic  and  we  ex¬ 
pect  about  half  to  be  partitioned  to  potential  energy  as  the  calcula¬ 
tion  proceeds  and  equilibrium  is  attained. 

The  computer  program  selects  initial  velocities  by  using  a  random 
number  generator  to  choose  the  initial  Ax,  Ay,  Az  for  each  particle  to 
be  between  -1  and  +1.  These  velocities  are  then  scaled  [see  Eq.  (109)] 
so  that 

N 

KE  -  \  m(ii  +  y\  +  *i>  -  2  (f  kT)  »  (149) 

i-1 

where  the  factor  2  is  included  according  to  the  discussion  above. 

The  initial  velocity  distribution  is  not  normal.  We  expect  that,  as 
the  calculation  proceeds  and  the  system  approaches  equilibrium,  the 
distribution  will  become  normal.  This  helps  determine  whether  or  not 
the  system  has  attained  equilibrium  and  will  be  discussed  in  Sec.  IV. D. 


-  ^  . 
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Note  also  that  the  temperature  is  not  specified  precisely  as  an  initial 
condition.  We  specify  the  initial  conditions  as  outlined  above,  and 
when  the  system  equilibrates,  we  calculate  the  system  temperature  from 
the  average  of  the  kinetic  energy  as  given  in  Eq.  (109). 

2.  hep  Initial  Positions.  We  specify  the  initial  particle  posi¬ 
tions  for  hexagonal  close-packed  crystals  by  placing  the  particles  at 
the  perfect  hep  lattice  sites.  Figure  13  shows  an  hep  lattice,  and  the 
lattice  vectors  a,  b,  and  c  in  the  Cartesian  coordinate  directions. 

With  these  lattice  vectors  there  are  four  particles  per  unit  cell, 
placed  at  positions  given  in  Table  I. 

These  unit  cells  are  repeated  throughout  the  calculational  volume 
to  yield  an  hep  structure  with  one  set  of  close-packed  planes  normal  to 
the  z-axis.  For  a  perfect  hep  structure,  which  corresponds  to  a  close¬ 
packing  arrangement  of  spheres,  the  relationships  between  the  lattice 

2 

vector  magnitudes  and  the  volume  per  particle  are 


and 


b  *  a/3  ,  c  *  a/8/3  , 


0  /3  2  a3  . 

■  -7-  a  c  »  — 

0  A  /2 


(150) 


TABLE  I 

PARTICLE  POSITIONS  WITHIN  AN  hep  UNIT  CELL 
Particle  x  y  z 


0 

1/2  a 
1/2  a 


0 

1/2  b 
1/6  b 


0 

0 

1/2  c 


ir 

4 


4 


0 


2/3  b 


1/2  c 
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Fig.  13. 

The  hexagonal  close-packed  structure  with  lattice 
vectors  2,  b,  and  c  indicated.  The  four  particles  in 
a  unit  cell  are  numbered. 
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The  radial  distances  to  the  shells  of  particles  relative  to  the  origin 
placed  at  one  of  the  lattice  sites  and  the  number  of  particles  within  a 
shell  are  shown  for  a  perfect  hep  lattice  in  Fig.  14.  The  bcc  posi¬ 
tions  are  also  shown  on  this  figure,  along  with  a  plot  of  the  effective 
pair  potential. 

An  hep  crystal  produced  by  the  molecular  dynamics  program  is  shown 
in  Fig.  15.  In  this  figure  the  near  neighbors  within  each  close-packed 
plane  normal  to  the  z-axis  have  been  connected  by  lines  for  clarity. 

The  dashed  and  dotted  vertical  lines  indicate  the  relative  positions  of 
the  planes.  Figure  16  shows  two  such  planes  as  viewed  looking  down  the 
z-axis.  Here  we  have  noted  the  traditional  A,  B,  C  designations  for  the 
relative  positioning  of  the  planes.  For  an  hep  structure,  the  close- 
packed  planes  are  stacked  in  an  ABAB* • *  arrangement.  For  a  face-cen¬ 
tered  cubic  (fee)  structure  the  stacking  is  ABCABCABC • • • . 

3.  bcc  Initial  Positions.  The  lattice  positions  of  a  body-cen¬ 
tered  cubic  structure  with  an  af  lattice  constant  (cube-edge  dimension) 
are  determined  easily.  The  distance  between  a  particle  and  its  nearest 
neighbor  is  and 

Rx  -  Y  a’  .  (151) 

There  are  two  particles  per  unit  cell  so  that  the  volume  per  particle 
is 


The  position  vectors  $(N)  to  all  points  in  the  lattice  may  then  be 
2 

specified  by 


Fig.  14. 

Radial  positions  of  particle  shells  for  the  hep  1 

and  bcc  crystal  structures  appropriate  to  a  volume 

per  particle  of  256  aj*  The  number  in  parentheses  \ 

is  the  number  of  particles  in  each  shell.  The  j 

total  ion- ion  interaction  potential  of  Fig.  9  is 

also  plotted.  ij 

i 

] 

1 


i 

< 


Fig.  15. 

The  hexagonal  close-packed  structure  as  set  up  by  the 
molecular  dynamics  program.  Solid  lines  are  drawn 
between  nearest  neighbors  in  each  close-packed  plane 
normal  to  the  z-direction  for  clarity.  The  dashed 
lines  indicate  the  relative  positions  of  the  planes 
marked  by  A  and  B. 


Fig.  16. 

Close-packed  hep  planes  viewed  down  the  2-axis.  The 
dashed  lines  indicate  an  A-plane  In  Fig.  15.  The  solid 
lines  indicate  a  B-plane.  For  an  hep  structure,  planes 
which  occupy  the  C  positions  do  not  occur. 
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R1 

R(N)  -—(».*+  N,?  +  N,£)  ,  (153) 

/3  1  L  3 

where  N^,  N ^ are  integers  and  are  constrained  to  be  either  all  even 
or  all  odd.  This  constraint  is  given  equivalently  by  requiring  that 
(Ni  +  N^).  (N^  +  N^),  and  (N^  +  N^)  be  all  even.  The  distance  to  each 
particle  is  given  by 

.  R,  ,  2  9  1/2 

|R(N) |  -  -±  (NT  +  N,  +  NT)  .  (154) 

/3  1  2  3 

Using  these  relations  we  easily  generate  the  points  in  a  perfect 
bcc  lattice  with  a  computer  program  and  we  have  used  this  method  to 
study  perfect  lattice  calculations  of  the  total  crystal  potential,  as 
described  in  Sec.  IV. C. 

However,  for  the  molecular  dynamics  calculations  we  -Hoose  to  use 

a  different  but  equivalent  method.  We  want  to  create  a  bcc  structure 

that  is  oriented  to  resemble  as  closely  as  possible  the  hep  structure 

that  we  are  studying.  We  will  create  a  bcc  system  with  (110)  close- 

packed  planes  normal  to  the  z-direction  because  our  hep  close-packed 

* 

planes  are  normal  to  the  z-direction. 

This  may  be  accomplished  by  setting  up  a  face-centered  tetragonal 
(fet)  lattice  with  lattice  vector  magnitudes  given  by 

a«c»/2b«/2»'  ,  (155) 

where  a*  is  the  desired  bcc  lattice  constant.  This  structure  is  shown 
schematically  in  Fig.  17.  The  dashed  lines  outline  the  cubic  box. 

*Based  on  a  procedure  by  Brad  Lee  Holian,  Los  Alamos  National  Labora¬ 
tory. 
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Fig.  17. 

A  face-centered  tetragonal  structure  with  a  -  c  - 
/l  b  appropriate  to  produce  a  body-centered  cubic 
(bcc)  structure  with  (110)  planes  normal  to  the 
z-direction.  The  four  points  in  a  unit  cell  are 
numbered.  The  dashed  lines  indicate  the  bcc  struc¬ 
ture. 
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The  bcc  (110)  planes  are  the  faces  of  the  fct  structure  normal  to  the 
z-axis . 

This  crystal  structure  has  four  particles  per  unit  cell  with 
positions  given  by  Table  II  and  indicated  in  Fig.  17. 

Figure  18  shows  the  bcc  lattice  as  set  up  by  the  molecular  dy¬ 
namics  program.  In  this  figure  a  body-centered  cube  is  outlined  and 
lines  are  drawn  through  the  cube  diagonal  and  body-centered  particle  to 
indicate  the  (110)  plane.  Figure  19  is  a  drawing  of  the  same  lattice 
with  lines  drawn  connecting  the  nearest  neighbors  within  a  close-packed 
plane.  The  numbers  1  and  2  designate  relative  positions  of  one  plane 
with  respect  to  another.  Figure  20  shows  a  "l"  plane  and  a  "2"  plane 
as  viewed  looking  down  the  z-axis.  These  last  two  figures  may  be  com¬ 
pared  with  Figs.  15  and  16  for  the  hep  lattice.  We  note  here,  in  pass¬ 
ing,  that  a  slight  compression  of  the  planes  in  Fig.  20  in  the  y-direc- 
tion  to  form  hexagonal  planes  and  a  relative  shift  between  planes  1  and 
2  in  the  x-direction  will  create  a  hexagonal  close-packed  structure  (or 
fee  structure,  depending  on  the  stacking).  The  radial  distances  to  the 
shells  and  the  number  of  particles  in  each  shell  for  a  bcc  lattice  are 
shown  relative  to  hep  structure  and  the  effective  pair  potential  in  Fig . 
14. 


TABLE  II 


PARTICLE  POSITIONS  WITHIN  THE  fct  UNIT  CELL 

Particle  x  y  z 

1 

0 

0 

0 

2 

1/2  a 

1/2  b 

0 

3 

1/2  a 

0 

1/2  c 

4 


0 


/  2  b 


1/2  c 
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Fig.  18. 

The  body- centered  cubic  structure  as  set  up  by  the 
molecular  dynamics  program.  The  lines  outline  the 
basic  cube  and  indicate  the  close-packed  plane 
through  the  cube  diagonal  and  body-centered  particle* 


Fig.  19. 

The  same  body-centered  cubic  structure  as  in  Fig.  18 
but  with  lines  drawn  between  the  nearest  neighbors 
within  a  close-packed  plane  for  comparison  with  the 
hep  structure  of  Fig.  15. 
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Fig,  20. 

Close-packed  planes  of  a  bcc  lattice  viewed  down  the 
z-axis.  The  dashed  lines  and  open  circles  Indicate 
planes  marked  1  In  Fig.  19.  The  solid  lines  and  solid 
circles  indicate  planes  marked  2. 
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C.  Determination  of  Run  Parameters 

In  this  section  we  determine  run  parameters  for  use  as  input  to  the 
molecular  dynamic  calculations.  These  parameters  are  the  time  step  At  of 
Eqs.  (90)-(92),  the  maximum  range  of  the  potential  RMAX  (see  Fig.  1),  and 
the  system  size.  There  is  no  standard  procedure  for  choosing  these 
parameters  and  generally  each  must  be  investigated  to  minimize  the  ef¬ 
fects  on  the  calculations.  We  discuss  how  we  have  determined  each  of 
these  parameters  for  use  in  our  calculations. 

1.  Determination  of  the  Time  Step.  The  time  step  must  be  kept 

3 

small  enough  that  errors  of  order  (At)  inherent  in  the  central  differ¬ 
ence  scheme  described  in  Sec.  II. B  are  negligible.  Another  way  of  think¬ 
ing  is  that  we  must  not  let  the  particles  in  the  system  move  very  far 
before  stopping  and  recalculating  the  forces  on  them.  Also,  we  must  not 
have  At  so  small  that  the  calculations  require  an  inordinate  amount  of 
computer  time. 

To  evaluate  At,  we  investigate  the  environment  that  a  single  par¬ 
ticle  in  our  system  experiences.  We  do  this  by  performing  a  "cell  model" 
calculation  where  one  particle  is  allowed  to  move  in  the  force  field  of 
all  the  others  which  are  constrained  to  their  perfect  lattice  positions. 
Doing  this,  we  can  map  the  potential  surface  for  a  particle  in  this  sys¬ 
tem. 

We  set  up  a  bcc  lattice  with  near  neighbor  distance  of  6.93  a^, 

3 

which  is  appropriate  to  an  atomic  volume  for  sodium  of  256  a^.  A  plot 
of  the  potential  surface  for  a  particle  moving  in  the  (001)  (or  xy) 
plane  in  the  +x  and  +y  directions  is  shown  in  Fig.  21.  This  ceil  model 
potential  well  is  quite  harmonic  and  symmetric,  as  shown  in  Fig.  22, 
where  the  potential  of  a  particle  moving  along  three  crystal  directions 


is  plotted  versus  the  square  of  its  distance  from  its  perfect  lattice 
site.  The  directions  were  chosen  to  indicate  the  most  and  least  dras¬ 
tic  paths  for  the  particle.  The  three  directions  shown  are  (1)  the 
[111]  direction  which  is  toward  the  body  center  (the  nearest  neighbor), 
(2)  the  [110]  direction  which  is  across  a  cube  face  diagonal,  and  (3) 

the  [100]  direction  which  is  down  a  cube  side. 

2 

The  curves  are  linear  in  r  out  to  at  least  an  atomic  unit  and  we, 
therefore,  represent  the  potential  in  this  region  by  a  harmonic  poten¬ 
tial  of  the  form 


$(r)  *  j  kr2  4*  C  . 


We  calculate  the  slope  of  the  line  for  the  [111]  direction  (the 
largest  slope  of  the  three)  and  find  that 


k  ^  0.0072  Ry/a^  , 


which  allows  us  to  calculate  the  period  of  a  harmonic  oscillation  in 
this  potential  as 


T  »  27tJ|  —  74  T0 


where  T0  is  the  time  unit  evaluated  in  Sec.  IV, A.  We  now  use  a  conven¬ 
tion  which  is  based  on  experience  and  says  that  a  conservative  time 
step  estimate  should  allow  a  particle  to  move  l/60th  of  its  period  per 
time  step.  Using  this  estimate  would  yield  a  At  of  about  1.0  T0.  We 
will  use  this  number  for  our  calculations. 

2.  Determination  of  RMAX.  We  choose  the  range  of  the  potential 
RMAX  (see  Fig.  1)  with  two  thoughts  in  mind.  We  want  the  effect  of  all 
particles  farther  away  from  a  given  particle  than  RMAX  to  be  negligible 
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and  we  want  to  minimize  the  effect  of  particles  moving  from  outside  to 
inside  the  RMAX  range.  Noting  the  small  magnitude  of  the  oscillations 
in  the  potential  and  force  functions  of  Figs.  9  and  10,  we  intuitively 
feel  that  an  RMAX  greater  than  16  a^  would  be  suitable.  To  get  a 
better  indication  of  the  effect  of  RMAX  on  the  calculations,  we  will 
calculate  the  perfect  crystal  total  potential  per  particle  for  the  be c 
lattice  for  different  values  of  RMAX. 

The  total  crystal  potential  per  particle  is  given  by  Eq.  CL06)  as 


,  N 

'■5U*V  • 

1-1  j-1 


(158) 


where 


r 


ij 


(159) 


To  calculate  this  for  a  perfect  crystal,  we  arbitrarily  place  one  par¬ 
ticle  at  the  origin  (particle  i)  so  that 


i-i 1  j-i 


*(rij} 


(160) 


The  sum  over  i  is  N  because  the  result  for  all  particles  in  a  perfect 
lattice  is  identical.  We  take  the  sum  on  j  to  be  over  all  the  other 
particles.  The  result  is 


$ 


N 

iZ*(r 

j-2 


(161) 


To  calculate  the  crystal  potential,  it  is  inappropriate  to  extend 
the  discrete  sums  to  infinity.  Additionally,  the  values  for  the 


molecular  dynamics  system  will  be  summed  discretely  using  Eq.  (106)  and 
will  be  cut  off  at  RMAX,  where  RMAX  must  be  kept  to  a  reasonable  value 
so  that  the  length  of  the  calculations  does  not  become  prohibitive. 

We  will,  therefore,  separate  the  discrete  sum  into  a  discrete  part 
(subscript  d)  out  to  RMAX  and  a  continuous  part  (subscript  c)  from  RMAX 
to  infinity,  so  that 


,  (162) 

a  c 

where  the  discrete  part  is  given  by 

r4  .<RMAX 
iJ 

♦d-7  E  ♦<'«>  •  'l63> 

j-2 

We  calculate  the  continuous  part  by  assuming  that  the  density 
approaches  a  uniform  distribution  at  large  r.  For  a  uniform  density 
with  one  particle  assigned  to  each  volume  per  particle,  ftg,  the  system 
potential  is  given  by 

/ao 

dr  r24>(r)  .  (164) 

-  0 

The  continuous  part  in  Eq.  (162)  is  given  by 


$ 

c 


-  f° 

Qo  1 


dr 


r20(r) 


RMAX 


(165) 


The  total  number  of  particles  within  a  sphere  of  radius  r  is 
plotted  in  Fig.  23  for  the  hep  and  bcc  crystal  structures  and  the  uni- 

3 

form  density  distribution  as  given  by  (47r/3^)r  .  Also  shown  in  this 
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Fig.  23. 

Number  of  particles  within  a  given  radius  vs  radius  for  the  hep 
and  bcc  crystal  structures  and  for  a  uniform  density  distribu¬ 
tion.  The  vertical  lines  on  the  r-axis  indicate  the  positions 
of  the  crystal  structure  shells  and  the  height  of  these  lines 
indicates  the  number  of  particles  in  each  shell. 


figure  are  the  positions  and  number  of  particles  in  each  shell  of  the 
crystal  structures. 

To  calculate  the  continuous  contributions  to  the  potential,  we  will 
use  the  large  r  limit  for  the  potential.  As  discussed  in  Sec.  III.B.5, 
in  this  limit  the  potential  goes  as 


cos  2k^r 
(2kfr)3 


We  calculate  this  term  to  be 


I  large  r 


*A(r) 


-!3  TTZ2*4f 


cos  2k^r 
(2kfr)3 


where  is  the  electron  kinetic  energy  at  the  Fermi  sphere  and 

is  the  magnitude  of  w  / e (q)  evaluated  at  q  *  2k,..  For  an  atomic  volume 

i5q  t 

3  -1 

of  256  a^,  which  yields  a  Fermi  wavenumber  of  0.4872  a^  ,  the  value  for 

w^^  is  0.0067.  We  call  $A(r)  the  asymptotic  form  of  the  potential  func¬ 
tion.  This  asymptotic  form  is  plotted  in  Fig.  24.  The  dashed  line  in 
this  figure  is  the  actual  potential  function  [0 (r )  of  Eq.  (77)]. 

Using  the  asymptotic  form,  we  calculate  the  continuous  part  of  the 
crystal  potential  to  be 


■••(*)( 


-Ci(2kf  •  RMAX) I 


where  Ci  is  the  cosine  integral  as  defined  in  Ref.  22. 

We  calculate  the  total  crystal  potential  using  Eqs.  (162),  (163),  and 
(167).  The  results  are  plotted  in  Fig.  25.  Large  jumps  are  noted  in  the 
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value  as  shells  of  atoms  in  the  structure  are  included  in  regions  where 
the  magnitude  of  the  potential  function  is  appreciable.  The  results 
shown  in  this  figure  indicate  that  a  choice  of  RMAX  greater  than  16  a^ 
is  suitable  because  the  magnitude  of  the  potential  is  small  and  the 
effect  of  particles  outside  RMAX  can  be  suitably  accounted  for  by  the 
continuous  contribution  $  .  There  is  one  more  point  to  consider  before 
investigating  the  magnitude  of  this  contribution. 

As  mentioned  at  the  beginning  of  this  discussion,  RMAX  must  also 
be  chosen  to  minimize  the  effect  of  particles  moving  from  outside  to 
inside  the  RMAX  range.  This  is  done  by  choosing  RMAX  to  coincide  with 
one  of  the  zeros  of  the  force  function.  In  this  way,  a  particle  sitting 
at  RMAX  would  have  seen  zero  force  whether  or  not  the  potential  was  cut 
off  there. 

With  the  above  discussion  in  mind,  we  have  somewhat  arbitrarily 
chosen  RMAX  to  be  at  the  zero  of  the  force  function  after  one  positive 
hump  and  one  negative  hump  (see  Fig.  10).  This  occurs  at  RMAX  *  21.65 

3 

a^  for  sodium  with  a  volume  per  particle  of  256  a^.  At  this  RMAX  we 
calculate  the  following  values  for  $,  and  their  ratio. 

$  -1.154  x  10-2 

$  1.1  x  10~5 

$  /#  9-5  x  10~4 

c 

The  continuous  contribution  to  the  potential  is  less  than  1%  of  its 
value.  Because  this  value  is  independent  of  the  details  of  the  molec¬ 
ular  dynamics  system  structure  and  is  much  smaller  than  the  absolute 
accuracy  of  our  calculations  (see  discussion  at  the  end  of  Sec.  V),  we 
will  neglect  it  for  this  study. 
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Because  the  potential  and  force  functions  depend  on  the  atomic  vol-  , 

ume,  the  appropriate  value  for  RMAX  will  also  vary.  Figure  26  is  a  plot 
of  the  RMAX  value  appropriate  for  each  atomic  volume  in  our  range  of 

interest.  A  linear  fit  of  the  points  of  this  plot  yields  > 

RMAX  -  0.008063  8Q  +  19.59  .  (168) 

We  use  this  value  of  RMAX  for  our  molecular  dynamics  calculations. 

3.  Determination  of  System  Size.  The  system  size  is  determined  in 
the  molecular  dynamics  program  by  specifying  the  number  of  unit  cells  to 
be  stacked  in  a  given  direction.  L^,  L^,  and  L ^  are  integers  that  spec¬ 
ify  the  number  of  unit  cells  in  each  Cartesian  direction,  respectively. 

If  a,  b,  and  c  are  the  magnitudes  of  the  lattice  vectors  in  these  direc¬ 
tions,  then  the  lengths  of  the  sides  of  the  computational  box  are 

xL  ®  L  -  a  , 
x 

yL  *  Ly  *  b  ,  and  (169) 

zL  -  L  •  c 
z 

We  are  constrained  by  the  minimum  image  convention,  discussed  in 
Sec.  II. B. 3,  so  that  each  box  length  must  be  at  least  twice  the  magni¬ 
tude  of  the  range  of  the  potential,  RMAX.  We  mentioned  in  Sec.  IV. B. 2 
the  possibility  of  the  hexagonal  close-packed  planes  normal  to  the  z- 
axis  being  stacked  either  in  the  ABAB***  (for  hep)  or  ABCABC #  *  *  (for 
fee)  arrangement.  Therefore,  we  would  like  to  have  a  multiple-of-six 
number  of  close-packed  planes  in  the  z-direction  so  that  neither  of 
these  possibilities  is  prohibited  by  the  periodic  boundary  conditions. 

With  these  constraints  in  mind,  we  choose  for  the  hep  structure 

L-7,L-4,L-6,  (170) 

x  y  z 


—  -*  *• 
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and  for  the  bcc  structure  (in  reality  a  face-centered  tetragonal  as 
discussed  in  Sec.  IV. B. 3) 

Lx  -  4  ,  Ly  -  7  ,  Lz  »  6  .  (171) 

Because  each  of  these  has  4  atoms  per  unit  cell,  they  contain  N  *  672 
particles  and  have  12  close-packed  planes  normal  to  the  z-axis.  Figures 
27  and  28  are  schematics  of  these  structures  in  their  initial,  perfect 
crystal  conf iguration. 

To  investigate  the  system  size  dependence,  we  have  calculated  the 
temperature  for  bcc  sodium  for  system  sizes  N  -  672,  864,  1372,  and 
2048  at  input  temperatures  of  50  and  300  K.  The  results  of  these  calcu¬ 
lations  are  shown  in  Fig.  29.  It  is  obvious  from  these  results  that  no 
definitive  statement  may  be  made  regarding  the  N  dependence  other  than 
it  is  small  and  within  +1 %  for  these  large  systems.  Therefore,  we  per¬ 
form  our  calculations  using  the  672-particle  systems  described  above. 

D.  Example  Calculations 

To  illustrate  the  molecular  dynamics  technique  and  the  approach  to 
equilibrium,  we  describe  two  calculations  in  detail.  The  calculations 
are  of  the  bcc  structure  described  in  Sec.  IV. B. 3,  performed  at  an  atomic 

3 

volume  of  256  a^  and  with  input  temperatures  of  50  K  and  300  K,  respec¬ 
tively. 

These  calculations  were  run  for  300  cycles,  which  corresponds  to  a 

-12 

time  of  150  time  units  or  1.05  x  10  seconds.  The  time  histories  of 
the  system  energies  per  particle  are  shown  in  Fig.  30.  The  solid  lines 
are  the  total  (potential  plus  kinetic)  energy,  the  dashed  lines  are  the 
potential  energies,  and  the  chain-dashed  lines  are  the  kinetic  energies. 


! 
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Fig.  27. 

a.  The  body-centered  cubic  structure  for  a  672-particle  system  as 
set  up  by  the  molecular  dynamics  program. 

b.  The  same  structure  but  with  nearest  neighbors  within  a  close- 
packed  plane  connected  by  lines. 
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It  is  obvious  from  these  figures  how  the  energies  equilibrate  from  the 
initial  conditions  where  the  particles  are  placed  at  perfect  lattice 
sites  and  assigned  initial  velocities  as  discussed  in  Sec.  IV.B.l.  As 
the  particles  move  from  this  initial,  unphysical  condition,  part  of  the 
energy  is  partitioned  to  potential  energy  and  the  system  equilibrates 
with  the  total  energy  remaining  constant. 

Once  the  system  relaxes  from  the  initial  conditions  and  begins  to 
oscillate  about  the  equilibrium  values  (this  happens  at  about  45  time 
units  for  these  calculations),  we  begin  our  time  average  and  average  to 
the  end  of  the  calculation. 

The  average  of  the  kinetic  energies  yields  the  temperature  by  way 
of  Eq.  (109).  The  average  total  energy  is  the  structure-dependent  en¬ 
ergy  of  the  system,  Eg,  which  will  be  added  to  the  volume-dependent 
energy  contributions,  as  discussed  in  Sec.  V,  to  yield  an  equation-of- 
state  point,  ETQT  (^q,T),  In  volume  and  temperature  space.  Table  III 
gives  the  Eg  and  T  values  obtained  from  these  calculations. 

An  indication  of  the  thermal  motion  of  the  particles  in  a  system 
may  be  obtained  by  calculating  their  atomic  distribution.  It  is  given 
by 

4irr2  d(r)  ,  (172) 


TABLE  III 

RESULTS  OF  THE  EXAMPLE  CALCULATIONS  FOR  bcc  SODIUM 
WITH  INITIAL  CONDITIONS  OF  50  K  and  300  K 


w 

T  00 

E 

s 

(Ry) 

256 

50.17  +  0.2 

-0.0105949 

+  1.  x  lo"7 

256 

293.38  +  0.6 

-0.0058439 

+  3.  x  10“7 
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where  d(r)  is  the  actual  particle  density  at  a  radius  r  from  a  reference 
particle.  We  calculate  this  distribution  for  each  particle  in  the  sys¬ 
tem  and  divide  the  total  by  the  number  of  particles  to  obtain  the  aver¬ 
age  atomic  distribution  for  the  system.  Figures  31a  and  31b  show  the 
atomic  distributions  for  the  systems  described  above.  The  solid  lines 
are  the  distributions  at  z  ^iven  instant.  The  dashed  lines  are  the 
distributions  averaged  ovsr  100  time  units,  which  indicate  the  perfect 
lattice  positions  about  which  the  particles  are  oscillating  (compare  with 
Fig.  14).  At  300  K  it:  is  difficult  to  distinguish  the  structure  of  the 
lattice. 

A  more  qualitative  indication  of  the  thermal  motion  of  the  par¬ 
ticles  is  seen  from  Figs.  32a  and  32b  which  are  plots  of  the  positions 
of  the  particles  at  a  given  instant  seen  looking  down  the  x-axis.  In 
the  perfect  lattice  position,  each  y,z  lattice  coordinate  would  show 
only  one  point.  At  300  K  there  are  indications  that  some  particles  may 
have  moved  out  of  their  perfect  lattice  positions.  We  will  investigate 
this  more  later  when  we  discuss  melting. 

As  mentioned  in  Sec.  II. B. 4,  an  indication  of  system  equilibration 
is  that  the  kinetic  temperatures  be  isotropic.  The  kinetic  temperatures 

are  T  ,  T  ,  and  T  ,  and  are  calculated  from  the  velocity  distributions 

x  y  z  7 

in  the  x-,  y-,  and  z-direc tions ,  respectively  [see  Eq.  (110)].  Figures 
33a  and  33b  are  time  history  plots  of  the  system  temperature  (solid  line) 
and  the  three  kinetic  temperatures  for  the  example  calculations.  The 
temperature  is  seen  to  be  isotropic  within  the  fluctuations  of  the  sys¬ 
tem  temperature. 

For  a  final  topic  regarding  these  two  example  calculations,  we 
discuss  a  measure  of  how  normal  the  velocity  distribution  is.  Such  a 
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Fig.  33. 

The  temperature  (solid  line)  and  kinetic  tempera¬ 
tures  (Tx,Ty,T2) (dotted,  chain-dotted,  dashed  lines) 
for  bcc  sodium  at  atomic  volume  of  256  ajj  and  tem¬ 
peratures  of  a)  50  K  and  b)  300  K. 


measure  is  the  kurtosis  [see  discussion  after  Eq.  (110)].  The  nth-order 
central  moment  of  the  N-particle  distribution  is  given  by 
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IE 


(173) 


where  x^  is  the  value  of  x  for  the  ith  particle  and  x  is  the  mean.  The 

kurtosis  measures  the  "degree  of  peaking"  of  the  distribution  and  may  be 

19 

defined  in  two  different  ways.  The  first  kurtosis,  denoted  by  3,  is 


6  -  —  •  (174) 

^2  is  the  second-order  central  moment  and  is  the  square  of  the  standard 
deviation,  a  *  (kT/m)  .  The  second  definition  of  kurtosis  we  denote  by 

c.23 

C  -  u4  -  3u2  .  (175) 

We  note  that 


6-3+-^-=3+ - - — -  .  (176) 

a  (kT/m) 

8  is  equal  to  three  for  a  normal  distribution  and  this  value  for  6  is 

19 

used  as  a  standard  to  indicate  how  normal  a  distribution  is.  C  is 
zero  for  a  normal  distribution. 

It  is  not  possible  to  precisely  predict  the  behavior  of  the  kur¬ 
tosis.  However,  we  may  derive  an  expected  dependence  from  a  simple 
model  and  then  compare  the  calculated  kurtosis  with  it  to  see  if  the  re¬ 
sults  are  consistent.  As  we  discussed  in  Sec.  II. B,  the  molecular  dy¬ 
namics  system  forms  a  microcanonical  ensemble  with  the  additional 
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constraint  that  the  system  momentum  is  conserved.  Holian  has  calculat¬ 
ed  the  kurtosis  for  a  microcanonical  ensemble  of  a  one-dimensional  chain 
of  hard  rods.  His  calculation  is  presented  in  Appendix  B.  For  large 
N,  he  predicts  that  the  kurtosis,  C,  is 


2 

so  that  we  expect  C  to  vary  as  T  /N.  Figure  34  is  a  time  history  plot 

2 

of  the  measured  kurtosis  C  times  the  factor  N/(kT/m)  for  the  50-K  and 

300-K  example  calculations.  We  see  that  C  rapidly  increases  from  its 

large,  negative,  initial  value  as  the  system  equilibrates.  The  value 
2 

for  CN/ (kT/m)  equilibrates,  roughly,  to  between  -100  and  0  for  both 
cases . 

With  these  numbers,  we  calculate  g. 


S 


CN/  (kT/tn) 2  13. 0 

J  N  12-85 


(178) 


The  measured  values  of  C  are  consistent  with  the  results  based  on  the 
simple  model  of  Appendix  B.  The  value  for  6  of  2.85  (close  to  3)  indi¬ 
cates  that  the  distribution  may  be  adequately  represented  by  a  normal 
19 

distribution. 

In  this  section  we  have  discussed  two  example  molecular  dynamics 
calculations.  We  have  found  that  the  system  we  are  using  is  adequately 
represented  by  a  normal  distribution  of  velocities  and  is  reasonably 
isotropic.  To  determine  the  equation-of-state  points  that  will  become 
our  data  for  studying  the  equation  of  state  of  sodium,  we  will  follow 

*Brad  Lee  Holian,  Los  Alamos  National  Laboratory. 
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the  procedure  outlined  above.  After  allowing  the  system  to  relax  from 
the  initial  conditions,  we  take  a  time  average  of  the  calculated  values 
to  yield  an  average  system  temperature  and  average  structure-dependent 
energy.  This  structure-dependent  energy  is  added  to  the  volume-depend¬ 
ent  terms,  as  described  in  Sec.  V,  to  yield  an  equation-of-state  point. 
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V.  CALCULATION  OF  THE  TOTAL  SYSTEM  ENERGY 

We  now  determine  the  expression  for  the  total  system  energy.  In 
developing  the  effective  ion- ion  interaction  potential  of  Sec.  II. A, 
we  neglected  the  terms  that  were  dependent  only  on  the  volume  of  the 
system  and  not  on  the  detailed  ion  arrangement.  We  will  here  evaluate 
these  terms  in  a  consistent  manner  so  that  they  may  be  added  to  the  en¬ 
ergy  calculated  by  the  molecular  dynamics  program  to  yield  the  total 


system  energy.  We  will  call  these  terms  'Volume-dependent  terms"  (E^) 


and  the  molecular  dynamics  energy  the  "structure-dependent  terms"  (E  ). 

s 


The  total  system  energy,  E^^,  is  given  by 


TOT 


T.  +  V.  +  E  *  Tt  +  V 
I  I  e  I 


(179) 


where  T^  is  the  kinetic  energy  of  the  ions,  is  the  ion- ion  potential. 
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and  Is  the  total  conduction  electron  energy.  We  write  V^.  as 


the  sum  of  the  coulomb  and  exchange  repulsion  terms  of  Eqs .  (75)  and 
(76). 


R  2L  r„  2  ^  E 


VT  “  + 

I  ES 


-Y„r 


B  ij 


(180) 


ij  « 


ij 


The  sums  are  over  all  ion  positions,  where  r^  is  defined  by 


rij  “  ,rj  -  v 


(181) 


and  the  factors  of  1/2  in  Eq.  (180)  take  into  account  double  counting. 

To  keep  track  of  the  volume-dependent  terms  in  the  total  conduction 


electron  energy,  E^,  we  evaluate  the  electron  energy  terms  of  Eq.  (13) 


in  a  slightly  different,  but  entirely  equivalent  manner.  We  rely  on 
the  local  pseudo potential  approximation  to  write  [see  Eqs.  (34),  (56), 
and  (57)] 


■jus- 

A 


4 

Hk 

3 


% 


i 


<icjw|jc  +  q>  »  W(r)e  ^  r  d^r  ■  W^  »  ^qWq 


We  use  the  subscript  0  to  denote  q  *  0  terms  and  we,  therefore,  write 
Eq.  (13)  for  the  electron  energy  in  the  kth  state  as 


%  •  S  +  “o -Z'  V-q(tk  -  W' 


where  the  prime  on  the  sum  indicates  that  the  q  *  0  term  is  omitted. 

We  have  constructed  the  local  pseudopotential,  which  consists  of  bare- 
ion,  screening,  and  exchange  and  correlation  terms  given  by  Eq.  (39)  as 


W-W_+W+W«W-+W  , 

B  s  x  B  sx  9 


where,  by  Eq.  (71) 


w  «  w  +  w 

B  z  c 


W^  is  the  coulomb  part  of  the  bar e-ion  potential  and  W^  is  the  core  re¬ 
pulsion  part. 

To  calculate  the  total  conduction  electron  energy,  we  must  sum 
over  all  the  occupied  states  and  subtract  off  the  electron-electron 
coulomb  self-energy,  E  ,  and  the  exchange  and  correlation  energies, 
Exc*  which  are  double-counted  in  the  sum  E^. 


The  coulomb  self-energy  is  given  by  Eq.  (52)  as 


Eee  "  f  E  dq  Vq  "  f  VsO  +  2  E  dq  Vq  * 


where  here  we  include  the  volume-dependent  q  *  0  term. 


—  ■  *  . 


i 


no 


In  dealing  with  the  exchange  and  correlation  effects,  we  have  as- 

2 

sumed  that  they  may  be  approximated  by  a  local  one-electron  potential 

2 

(see  Sec.  II. A. 5).  Wallace  discusses  the  approximation  for  exchange 
that  depends  only  on  the  density  of  the  conduction  electrons.  The  total 
exchange  energy  of  the  conduction  electrons  in  their  ground  state  is 
assumed  to  be 


k  J 


(187) 


where  X(r)  is  the  operator  that  represents  exchange  and  correlation  and 
is  a  function  of  the  density  only.  He  uses  a  variational  calculation  to 
relate  the  exchange  potential  Wx  to  X  by 


3(dX) 

dd 


(188) 


where  d  is  the  conduction  electron  density.  He  evaluates  the  exchange 

2 

and  correlation  double-counting  correction  as 

E*c ■-mo(xo - v  +!£' v,  •  <i8,) 

q 

so  that  the  total  double-counting  correction  is 

E  +  E  -  §  drtW  „  -  nd.CX.  -  W  n)  +  §  Y' d  W  .  (190) 

ee  xc  2  0  sO  0  0  xO  2  q  sx-q 

q 

We  may  now  write  the  total  conduction  electron  energy  as 

Ee  ’I)  W  -  (E,  +  E*c>  •  (191) 

k 


The  sum  over  the  occupied  states  is 


where  we  have  used  the  expression  for  the  Fourier  component  of  the 
density  of  Eq.  (45)  and  d^  ■  z/ft^.  We  arrive  at  the  following  ex¬ 
pression  for  the  total  conduction  electron  energy: 


■?  Vk + 


NzMrt  -  dnM  _  +  fidn(Xft  -  W  A)  +  d  (W  -W  )  . 

0  2  0  sO  0  0  xO  2  q  -q  sx-q 


At  this  point  we  must  digress  for  a  moment  to  discuss  how  the 
total  adiabatic  potential,  V  ■  +  E^,  is  evaluated  for  a  crystal 
structure.  From  Eqs.  (179),  (180),  and  (193),  V  is  given  by 


tl  £,  +  NzW  -  ^  d.W 
Tc  k  0  2  OsO 


VX0-“x0>+§E'd1("-o-,W  • 


q  sx-q 


The  last  term  in  this  expression  is  evaluated  using  the  relations  of 
Sec.  II. A  to  be 


Sv’d  (W  -  W  )  -  nT'  S  S  F(q)  -  E. 

2  ^  qx  -q  sx-q  L*  q  -q  n  bs 


This  is  the  band  structure  energy  of  Eq.  (30)  with  F(q)  given  by  Eq. 
(79)  so  that  the  expression  includes  the  screening,  exchange,  and  cor¬ 
relation  effects  as  discussed  in  Sec.  II. A.  Note  that  the  first  sum  in 
Eq.  (194)  is  divergent  as  are  the  Wg^  term  [see  Eq.  (53)]  and  the 


t 
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part  of  Wq  [see  Eq.  (72)].  A  standard  evaluation  of  the  first  sum  is  to 
2 

use  the  Ewald  method  whereby  the  first  sum  in  Eq.  (194)  is  split  into 

two  sums,  one  in  q  space  and  one  in  real  space.  The  divergent  term  in 

2 

this  sum  cancels  with  the  Wg^  and  W^  terms  so  that  the  expression  for 
V  is  finite.  The  band  structure  energy  is  evaluated  readily  for  a  per¬ 
fect  crystal  because  the  structure  factors  are  delta  functions  about  the 
reciprocal  lattice  vectors,  Q,  and 


which  yields  a  finite  result.  Wallace  discussed  this  method  in  detail 
and  used  it  to  evaluate  the  total  adiabatic  potential  for  sodium  and 
potassium. ^ 

For  our  purpose  here,  however,  we  are  calculating  the  terms  for  the 
adiabatic  potential  of  Eq.  (194)  differently  and  we  must  ensure  that  we 
are  properly  accounting  for  all  the  terms.  Having  restructured  the  q- 
space  sum  for  the  band  structure  energy  of  Eq.  (195)  in  terms  of  a  real 
space  sum  of  the  effective  interaction,  Vj^(r),  over  arbitrary  ion 
positions^’ ^  [see  Eq.  (31)],  we  note  an  obvious  difficulty  since  the 


E'wv  <1,7> 

ij 

is  divergent  because  of  the  leading  1/r  term  at  large  r  [see  Eq.  (139)] 
in  VIND(r)- 

The  difficulty  arises  because  the  q-space  sum  of  Eq.  (195),  origi¬ 
nating  from  the  perturbation  calculation  [see  Eq.  (13)],  explicitly 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 
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leaves  out  the  divergent  q  *  0  term,  yet  we  have  reintroduced  it  by 
transforming  to  the  sum  in  r  space,  which  we  take  to  infinity.  To 
correct  for  this  discrepancy,  we  allow  the  q  »  0  term  to  be  part  of 
the  q-space  sum  when  we  restructure 


as 


where 


NE  V-qF(q) 

q 


(198) 


fL’Tn«><rij>  • 


(199) 


lj 


IND 


<r>  • 


(200) 


which  are  the  same  as  Eqs .  (31)  and  (32)  with  the  primes  missing  from 
the  q  sums.  E^,  as  defined  here,  does  not  have  a  finite  value.  Note 
that,  with  the  help  of  Eq.  (195), 


so  that  we  may  write  for  the  total  adiabatic  potential 


* 


i 


tM***  *  *  f 
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2  2 


fZ'  rf  4E'v"VlJ  *Zv,*n-!v 

ij  iJ  ij  k 


sO 


*  CVX0  -  “*0>  -  2  VW0  -  "s,0>  *  I E  V%  -  W,x-,>  •  <2°2> 


The  last  sum  is  now,  in  a  consistent  way,  equal  to 

lZ*vim)(rij5  +£f(<i)  ■ 
ij  q 

The  terms  may  be  combined  and  arranged  to  yield 


(203) 


«_1V 


2e2  'YBrij  .  . 

_L.-  +  V  ♦*»*« 


-  v'f- 

2N  Lm*  r 

ij  L 

+ *(r  **<>-¥)  4  z  «■»  • 


+  N  Z  £ 


kwk 


(204) 


This  is  the  expression  for  the  total  system  potential,  which  is  added  to 
the  kinetic  energy  of  the  ions  to  yield  the  total  system  energy  of  Eq. 
(179). 


etot 


T.  +  V  *  T.  +  VT  +  E 
I  lie 


(205) 


The  term  in  the  square  brackets  of  Eq.  (204)  is  the  effective  ion- ion 
interaction,  $( r),  of  Eq.  (77)  and  the  1/r  coulomb  term  is  canceled  at 
large  r  by  the  leading  -I/r  term  in  vIND(r)*  as  discussed  in  Sec.  Ill; 
therefore,  all  the  terms  in  Eq.  (204)  are  finite. 

The  ion  kinetic  energy  is  the  time-averaged  molecular  dynamics 
kinetic  energy,  as  discussed  in  Sec.  II.  The  sum  in  the  brackets  of 
Eq.  (204)  is  the  time-averaged  molecular  dynamics  potential  energy. 


I 
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The  sum  of  the  two  is  the  total  structure-dependent  system  energy,  Eg, 
and  is  calculated  by  the  molecular  dynamics  program  as  discussed  in 
Sec.  III.B.5.  The  remaining  terms  of  Eq.  (204)  are  the  volume-depend¬ 
ent  terms,  Ev>  which  must  be  calculated  separately.  We  will  now  evalu¬ 
ate  these  terms. 

As  in  Sec.  III.A.l,  we  find  it  convenient  to  write  the  volume-de¬ 
pendent  terms  as  functions  of  rg,  the  radius  of  the  sphere  whose  volume 
is  the  average  conduction  electron  volume, 


4  3 

“t  Trr 
3  s 


1  r 

We  have  evaluated  —  l  n,  £,  in  Sec.  II. A. 2  as  an  integral  over  the 
w  k  K  K 

Fermi  sphere  to  arrive  at  the  average  kinetic  energy  of  the  electrons 
times  the  valence 


nS  nk£k 


z  }  ef 


2.210  „ 

z  — 5 —  Ry 

T 

S 


Wq  is  given  by  Eq.  (126),  so  that 


Z0  m  zz  m  -z(  1.22772)  Rv 
2  3  f  2  y 

rs 


To  evaluate  —  W  we  use  Eq.  (61)  to  write 


f  wxo  "  ^  =hST  f  (q)  dq 
q-K)  q 


With  Eq.  (81) 


2(q2  +  £k2) 


d°’S  * 


we  see  Chat 


zW  2  2 

_jtO  „  -w  e,  ,  .(0.407258)  £  7-  Ry  .  (212) 

The  remaining  volume-dependent  term  is  —  £  F(q),  which  we  convert 

N  q 

to  an  integral,  so  that 


q  F(q)  dq  . 


This  integral  is  evaluated  by  numerical  integration  using  the  same  pro¬ 
cedure  by  which  we  evaluate  the  integral  of  Eq.  (78),  as  described  in 
Sec.  III.B.5. 

Xq  is  the  total  exchange  and  correlation  energy  per  electron  for  a 
uniform  electron  gas,2 


Xn  *  £  +  £ 

0  X  C 


ex  is  the  Hartree-Foch  exchange  energy  and  is  calculated  to  be  the 
standard  result 


0.916  „ 

ze  *  -z  -  Ry 

x  r 


We  approximate  the  correlation  energy  per  electron,  using  the 
Pines-Nozieres  formula. ^  They  interpolate  between  expressions  that  are 
valid  in  low-q  and  high-q  regions  in  the  same  sense  as  the  Hubbard  in¬ 
terpolation  formula  discussed  in  Sec.  II. A. 5.  The  Pines-Nozieres  formu¬ 
la  yields 

ze  -  z(-0 .115  +  0.031  In  r  )  .  (216) 

c  s 

All  of  these  volume-dependent  terms  are  evaluated  separately  and 

added  to  the  molecular  dynamics  results.  Values  of  the  individual  terms 

3 

evaluated  at  an  atomic  volume  of  256  are  given  in  Table  IV.  Table  V 
gives  the  total  volume-dependent  contribution  to  the  system  energy  for 
several  different  volumes.  This  table  also  gives  the  results  of  the 
molecular  dynamics  calculation  of  the  static  (T  *  0)  potential  energy 
per  ion  for  bcc  sodium.  The  sum  of  the  two  terms  is  given  in  the  third 
column  and  is  the  total  static  potential  of  the  system  at  the  volume 
indicated.  The  values  of  this  table  are  plotted  in  Fig.  35,  where  the 
delicate  interaction  between  the  structure-  and  volume-dependent  terms 
is  evident. 

The  volume-dependent  terms  are  not  determined  absolutely  beyond  the 
second  or  third  decimal  place.  However,  we  will  be  observing  differ¬ 
ences  between  structures  at  the  same  volume  that  are  two  or  three  orders 
of  magnitude  smaller.  Therefore,  we  will  calculate  the  values  to  the 
accuracy  of  Table  V  and  recognize  their  validity  when  we  have  been 
treating  the  volume-dependent  terms  the  same  for  each  of  the  structures 
but  not  necessarily  in  an  absolute  sense. 

We  now  have  all  the  information  necessary  to  calculate  the  total 
system  energy.  We  will  perform  a  molecular  dynamics  calculation  at  a 


V 


TABLE  IV 

CALCULATED  VALUES  OF  THE  VOLUME- DEPEND ENT  TERMS  OF  EQ.  (214) 
EVALUATED  AT  AN  ATOMIC  VOLUME  OF  256  a\,  r  -  3.9390 


Expression 


ISni, 


/2.210\ 

■fai 

S 

-ft*) 


-2(0.115  -  0.031  in  r  ) 


/1.22772\ 

-<—) 

s 

z  /o.4Q7258\ 
rs  / 


27T2  J0  IT*™  ^ 


Value 

(Ry) 

0.14244 


-0. 23255 


-0.07250 


-0.07913 


0.05703 


Total 


-0.27828 


-0.46299 


After  the  original  publication  of  this  dissertation,  a  program¬ 
ming  error  was  found  in  the  subroutine  which  calculates  this 
integral.  Correction  of  this  error  results  in  a  difference 
in  the  calculated  volume-dependent  terms,  Ey,  of  about  -0.02  Ry, 
causing  the  reported  values  to  agree  better  with  experiment. 

The  new  values  have  been  included  here,  where  necessary,  or 
errata  notes  specified. 
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TABLE  V 

AT  THE  VOLUME  INDICATED,  COLUMN  2  IS  THE  SUM  OF  THE  VOLUME-DEPENDENT 
TERMS  OF  EQ.  (204),  COLUMN  3  IS  THE  STATIC  (T  *  0)  POTENTIAL  PER  ION 
FOR  bcc  SODIUM,  AND  COLUMN  4  IS  THE  TOTAL  STATIC  SYSTEM  ENERGY 


E 

E 

E  +  E 

0 

V 

s 

v  s 

lail 

(Ry) 

(Ry) 

(Ry) 

232 

-0.466836 

-0.0070983 

-0.473934 

240 

-0.465522 

-0.0087709 

-0.474293 

250 

-0.463924 

-0.0105869 

-0.474511 

256 

-0.462989 

-0.0115468 

-0.474536 

260 

-0.462374 

-0.0121400 

-0.474514 

270 

-0.460868 

-0.0134687 

-0.474337 

280 

-0.459405 

-0.0146061 

-0.474011 

given  volume,  and  input  temperature.  Once  the  system  has  equili¬ 
brated,  as  discussed  in  Sec.  IV. D,  we  will  calculate  the  average  tem¬ 
perature,  T,  potential  energy  per  particle,  and  kinetic  energy  per  par¬ 
ticle.  To  the  potential  and  kinetic  energies,  we  will  add  the  volume- 
dependent  energy  to  obtain  the  total  system  enerty,  E^^CQ^jT) ,  as  an 
equation-of-state  point  in  volume  and  temperature  space. 


-4.4H0.fJ7  -0.4H  -0.474  -0.474  -0.473  -a 473  -a 471  -0.470  -0.401  -O 


121 


! ; 


0 


VI.  RESULTS 

A,  Equation-of-State  Points  for  Solid  Sodium 

In  this  section  we  present  the  results  of  calculations  of  the  total 
system  energy,  E^^,  f°r  solid  sodium  in  the  hep  and  bcc  phases. 

1.  Static  Crystal  Potential.  We  determined  the  static  (zero  tem¬ 
perature)  system  energy  for  sodium  by  calculating  the  crystal  potential 
due  to  the  ion- ion  interaction  potential,  cj>(r) ,  using  Eq.  (106)  and  add¬ 
ing  to  it  the  volume-dependent  terms  as  discussed  in  Sec.  V.  The  results 
of  this  calculation  for  bcc  sodium  were  presented  in  Table  V  and  are 
plotted  as  the  lower  curve  in  Fig.  35.  The  results  for  both  the  bcc  and 
hep  structures  are  plotted  in  Fig.  36  as  a  function  of  volume.  We  will 
call  this  energy  the  static  crystal  potential 

We  performed  a  least-squares  fit  to  the  calculated  points  in  Fig.  36 
to  a  cubic  polynomial, 

W  *  Pj_  +  P2^0  +  P3fi0  +  P4n0  '  (217) 

The  coefficients  of  this  polynomial  are  given  in  Table  VI. 


TABLE  VI 

COEFFICIENTS  OF  THE  CUBIC  FIT  OF  THE  CALCULATED  STATIC  CRYSTAL 
POTENTIAL  AS  A  FUNCTION  OF  ATOMIC  VOLUME 


Term 

p2  (Ry) 
p2  (Ry/ao> 

P3  (Ry/a^) 
p4  (Ry/a^) 


_ bcc _ 

-0.310428  +  0.011 
(-1.676211  +  0.13)  x  10' 
(5.573895  +  0.51)  x  10‘ 
(-5.985510  +  0.66)  x  10‘ 


_ hep _ 

-0.297263  +  0.035 
(-1.828394  ±  0.42)  x  10' 
(6.161591  +  1.6)  x  10"' 
(-6.744432  +  2.1)  x  10~ 
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Fig.  36. 

The  total  static  crystal  potential  for  hep  (circles)  and  bcc  (squares) 
sod  ium. 
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For  this  classical  system,  at  zero  temperature,  the  Helmholtz 
free  energy,  F,  is  equal  to  the  static  crystal  potential  and  we  calcu¬ 
late  the  zero  temperature  pressure  from 


and  the  zero-temperature  bulk  modulus  from 


B 


0 


(218) 


(219) 


We  calculated  the  zero-temperature  pressure  and  bulk  modulus  using  the 

fits  to  Eq.  (217).  The  results  are  plotted  in  Figs.  37a  and  b. 

The  measured  zero-temperature  and  pressure  static  crystal  bulk  modu- 

21 

lus  for  bcc  sodium  are  reported  by  Wallace  as 


$0  “  -0.46  Ry  , 

=  255.5  ,  (220) 

and  Bq  =  5.05  x  10  ^  Ry/aj^ 

These  are  the  values  he  used  to  fit  the  parameters  a  ,  6,  and  p  in  the 

B 

pseudopotential,  as  discussed  in  Sec.  III. A. 2. 

For  bcc  sodium  our  calculated  is  -0.475  Ry  (see  Fig.  36).  We 

3 

calculate  the  zero-pressure  volume  to  be  255.1  a^,  which  i*  0.2%  different 
from  the  observed  value,  and  the  bulk  modulus  to  be  5.08  x  10  Ry/a^, 
which  is  0.5%  different  from  the  observed  value. 

As  is  apparent  from  Fig.  36,  we  calculate  that  the  hep  structure  is 
the  stable  phase  for  sodium  at  zero  temperature.  We  calculate  an  energy 


V 


125 

-5  3 

difference  of  14  x  10  Ry  at  an  atomic  volume  of  256  a^.  The  experi- 

24  -5 

mental  value  reported  by  Straub  and  Wallace  is  3.15  x  10  Ry,  and 

their  calculated  value  is  6  x  10  ^  where  they  have  included  the  quantum 

mechanical  zero-point  energy  of  the  sodium  crystal  structures. 

With  <J>q  established,  we  proceed  in  the  next  section  to  calculate 

equation-of-state  points  at  temperatures  greater  than  zero. 

2.  Total  System  Energy.  In  this  section  we  present  the  results  of 

calculations  of  the  total  system  energy,  E^^,  for  bcc  and  hep  sodium  at 

finite  temperatures.  The  calculations  proceed  as  discussed  in  Secs.  IV 

and  V.  For  each  calculation  the  volume  is  predetermined  as  an  input 

parameter,  the  temperature  is  calculated  from  the  time  average  of  the 

system  kinetic  energy,  and  the  structure-dependent  energy  Eg  is  the  time- 

averaged  sum  of  the  potential  and  kinetic  energies.  We  add  Eg  to  the 

volume-dependent  term  Ev  to  arrive  at  the  total  system  energy  ETqT  *  E  * 

E  +  E  . 
s  v 

We  know  from  harmonic  theory  that  the  system  energy  per  particle 
must  vary  as 

E  *  +  3kT  ,  (221) 

as  T  approaches  zero,  k  is  Boltzman's  constant.  We,  therefore,  write 
the  energy  as 


E(-Qq#T)  -  +  3kT  +  f (fl0,T)  ,  (222) 

so  that  the  function  f(.Qg,T)  contains  all  contributions  that  are  not 
harmonic.  Our  calculations  thus  become  a  direct  measure  of  this  function. 
The  results  of  the  calculations  of  the  total  system  energy  at  an 

3 

atomic  volume  of  256  a^  are  shown  in  Fig.  38.  The  solid  straight  lines 
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Fig.  38. 

Total  system  energy  vs  temperature  for  hep  and  bcc  sodium 
at  an  atomic  volume  of  256  a2. 
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are  drawn  from  for  each  structure  at  a  slope  of  3k.  The  error  bars 
on  the  calculations  are  less  than  the  size  of  the  circles  and  squares. 
The  deviation  from  harmonic  behavior  as  temperature  increases  is  obvious 
from  the  figure. 

We  calculated  the  system  energy  E(Qq,T)  for  bcc  and  hep  sodium  at 
four  different  volumes — 232,  250,  256,  and  270  a^ — and  at  temperatures 
between  0  and  400  K.  The  results  of  these  calculations  are  tabulated  in 
Appendix  C. 

From  the  calculated  values  of  E(.Qq,T),  $q(v),  and  T,  we  calculate 
the  function 


f(n0.T)  -  E(QQ,T)  -  *q(£20)  -  3kT 


directly.  The  results  of  this  calculation  are  plotted  in  Fig.  39.  The 
error  bars  assigned  to  the  values  are  the  calculated  standard  deviations 
of  the  means  of  the  statistical  time  averages.  The  major  portion  of  the 
error  in  the  calculated  point  is  due  to  the  statistical  fluctuations  of 
the  temperature. 

2 

We  fit  the  data  at  each  atomic  volume  with  a  T  form,  which  is  all 
the  accuracy  of  the  calculation  justifies.  The  function  f (0.^,1)  then 
takes  the  form 


fw0,T>  -  cayr  . 


The  values  of  C(Q^)  determined  by  a  least-squares  fitting  procedure  are 

given  in  Table  VII.  These  values  may  be  compared  with  the  value  of  the 
2 

T  coefficient  to  the  free  energy  of  sodium  estimated  from  heat  capacity 

2  -10 

measurements  and  reported  by  Wallace  as  8  x  10 


f(fl0,T)(Ry)  XKT*  f(flo.T)(Ry) 


128 


39a 


Fig.  39. 

The  function  defined  by  Eq.  (223)  plotted  vs  tempera  tore 

for  (a)  bcc  and  (b)  hep  sodium. 
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TABLE  VII 

VALUES  OF  THE  COEFFICIENT  C  (P.q) 


n 

o 

bcc 

VJ 

hep 

(x  10“9  Ry/k2) 

(x  10~9  Ry/k2) 

232 

0.908  +  0.64 

1.301  +  0.23 

250 

1.012+0.52 

1.475  +  0.50 

256 

1.250  +  0.10 

1.563  +  0.25 

270 

1.350  +  0.51 

2.057  +  0.33 

These  values  are  plotted  in  Fig.  40  along  with  a  linear  least- 
squares  fit  to  the  data.  The  values  for  the  fit 

c<V  -  cx  +  c2f<0  (225) 

are  given  in  Table  VIII. 

The  calculated  values  of  f<fl^,T>  for  the  bcc  structure  above  350  K 
in  Fig.  39a  were  not  used  for  the  fits  described  above.  The  significant 
deviation  of  the  270  a*  point  indicates  that  something  other  than  anhar- 
monic  effects  are  taking  place.  We  checked  this  point  and  found  that 
some  particles  were  diffusing  from  their  lattice  sites  so  that  partial 
melting  had  taken  place.  We  will  discuss  melting  in  the  following  sec¬ 
tions  . 


TABLE  VIII 

VALUES  FOR  THE  COEFFICIENTS  (Ry/k2)  AND  C2  (Ry/k2a2)  OF  EQ.  (225) 


1 

2 


_ bcc _ _ 

(-1.958  +  3.4)  x  10' 
(1.225  +1.3)  x  10' 


_  hep _ 

(-3.247  +  5.7)  x  10' 
(1.923  +  2.2)  x  10' 
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In  this  section  we  have  reported  the  results  of  the  calculation  of 
the  total  system  energy  of  solid  sodium  in  the  bcc  and  hep  structures. 

It  is  given  by 

E(fi0,T)  -  $0(G0>  +  3kT  +  C(£20)T2  ,  (226) 

which  is  valid  within  the  errors  of  our  calculations  for  atomic  volumes 
3 

from  230  to  279  a^  and  temperatures  from  0  to  300  K. 

B.  The  Liquid  Phase 

We  can  cause  the  molecular  dynamic  system  to  melt  by  beginning  a 

calculation  at  a  temperature  high  enough  that  the  particles  have  enough 

velocity  to  move  from  their  crystal  lattice  positions  and  diffuse  through 

3 

the  system.  We  did  this  for  bcc  sodium  at  an  atomic  volume  256  a^  by 
choosing  an  input  temperature  of  700  K.  The  temperature  history  of  this 
calculation  is  shown  in  Fig.  41a.  If  the  system  were  harmonic  the  tem¬ 
perature  would  begin  to  oscillate  about  700  K  because  the  initial  energy 
would  be  partitioned  equally  between  the  potential  and  kinetic  energies 
(see  discussion  in  Sec.  IV.B.l).  Due  to  the  anharmonicity  of  the  system, 
the  initial  temperature  is  about  620  K  (see  Fig.  41a).  As  diffusion 
occurs,  the  potential  energy  of  the  system  is  increased  by  particles  mov¬ 
ing  out  of  their  low-potential  energy  lattice  positions.  This  occurs  at 
the  expense  of  the  kinetic  energy  because  the  total  energy  in  a  molecular 
dynamics  calculation  is  conserved.  The  potential,  kinetic,  and  total 
energies  are  shown  versus  time  in  Fig.  41b.  Thus  the  temperature  of  the 
system  is  lowered  and  the  equation-of -state  point  moves  horizontally  to 
the  left  in  an  E,T  plot,  as  Indicated  by  the  arrow  at  point  1  In  Fig.  42. 

We  then  use  this  calculation  as  a  starting  point  to  calculate  equa- 
tion-of-state  points  at  other  temperatures.  We  do  this  by  artificially 


multiplying  the  velocity  of  each  particle  by  a  factor  close  to  1.0  to 
slowly  decrease  or  increase  the  temperature  of  the  system.  Then  when 
the  system  has  reached  the  temperature  at  which  we  want  the  equation- 
of -state  point,  we  set  this  factor  equal  to  1.0  and  allow  the  system  to 
equilibrate  and  take  time  averages  just  as  we  did  for  the  equation-of- 
state  points  arrived  at  in  Sec.  VI. A. 

In  this  manner  we  investigated  the  equation  of  state  for  liquid 

3 

sodium  at  an  atomic  volume  of  256  a^.  The  results  are  the  square  points 
shown  in  Fig.  42.  The  circles  are  the  calculated  bcc  solid  points  which 
are  plotted  here  for  comparison.  The  calculations  were  performed  in  the 
order  indicated  by  the  numbers  on  the  points.  The  system  came  down  the 
solid  curve,  up  the  lower  dashed  curve,  down  the  upper  dashed  curve,  and 
back  up  the  lower  dashed  curve.  The  hysteresis  that  is  evident  in  this 
figure  Is  due  to  the  system  ’’freezing  out"  available  lower-energy  con¬ 
figurations  as  the  temperature  is  lowered. 

The  interesting  point  here  is  that  the  liquid  phase  defines  its 
equation-of-state  curve  back  to  low  temperatures  where  it  has  formed  a 
metastable,  glassy  system.  Within  the  confines  of  the  periodic  system 
and  for  the  short  times  (in  a  physical  sense)  of  a  molecular  dynamics 
calculation,  this  glassy  state  is  stable  enough  to  calculate  equation-of- 
state  points  as  we  did  for  solid  bcc  and  hep  sodium.  We  recognize,  of 
course,  that  at  low  temperatures  the  bcc  phase  of  sodium  is  also  a  metast¬ 
able  state  with  respect  to  the  hep  phase. 

This  glassy  state  is  most  unstable  at  low  temperatures.  We  allowed 
the  system  initially  at  5  K  to  equilibrate  for  1200  cycles.  During  this 
time  it  increased  its  temperature  to  9  K  by  adjusting  to  a  lower 
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potential  energy  configuration  and  moved  along  the  line  marked  ’’DRIFT1* 

in  Fig.  42,  Twelve  hundred  cycles  is  a  rather  long  molecular  dynamics 

-12 

calculation,  but  represents  a  very  short  real  time  (4  x  10  s) .  We 
will  calculate  the  equation  of  state  of  liquid  sodium  by  starting  at  a 
low- temperature,  glassy  state  as  the  initial  configuration,  raising  the 
temperature,  and  allowing  the  system  to  equilibrate,  just  as  we  did  for 
the  solid  systems. 

Because  the  glassy  system  drifts  noticeably  at  low  temperatures,  we 
must  estimate  the  static  potential  of  this  state.  We  do  this  by  calcu¬ 
lating  the  metastable  equilibrium  at  a  low,  but  finite,  temperature  (ap¬ 
proximately  30  K) .  At  this  temperature  we  expect  that  the  system  is  es¬ 
sentially  harmonic  in  that  the  system  energy  will  be  equally  partitioned 
between  kinetic  and  potential  (so  that  APE  *  KE  in  Fig.  43).  We  then 
estimate  the  static  potential  for  the  glassy  state  as 

*  PE  -  KE  .  (227) 

We  changed  the  volume  of  the  glassy  state  by  changing  the  periodic 
box  dimensions  and  scaling  the  particle  positions  in  the  same  ratio  so 
that  each  particle  retained  the  same  relative  position  within  the  box. 

We  then  estimated  the  static  potential  using  Eq.  (227).  The  results  are 
plotted  in  Fig.  44. 

We  fit  the  calculated  points  with  a  cubic  polynomial, 

W  -  PX  +  P2«0  +  p3^0  +  p4fi0  *  (228) 

where  the  coefficients  have  the  following  values. 
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Fig.  43. 

Time  histories  of  the  kinetic,  potential,  and  total 
energies  of  a  molecular  dynamics  calculation.  If  the 
system  is  harmonic,  APE  *  KE. 


Fig.  44. 

The  estimated  static  potential  for  the  glassy  state  of  sodium 
plotted  vs  atomic  volume. 
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Px  -  -0.304076  +  0.0036 
P2  -  (-1.712226  +  0.065)  x  l(f3 
P3  -  (5.660021  +0.26)  x  10_6 
p4  -  (-6.067466  +  0.33)  x  10_9 

We  calculated  the  finite  temperature  equation-of-state  points  for 

3 

the  liquid  state  of  sodium  at  atomic  volumes  232,  256,  and  270  a^.  The 
results  of  these  calculations  are  tabulated  in  Appendix  C  and  plotted  in 
Figs.  45a,  b,  and  c,  respectively.  The  equation-of-state  points  for  bcc 
sodium  at  these  volumes  are  included  in  these  figures.  The  circles  denote 
points  that  were  calculated  starting  with  the  bcc  configuration  and  the 
squares  denote  points  that  were  calculated  starting  at  the  glassy  state. 
The  circled  points  on  the  liquid  curve  have  melted  from  the  bcc  config¬ 
uration. 

We  analyzed  the  equation-of-state  points  for  the  glassy  state  in 
much  the  same  way  we  analyzed  the  solid  equation-of-state  points  in 
Sec.  VI. A.  As  before,  we  expect  the  temperature  dependence  of  the  energy 
to  be  of  the  form 


E«20,T)  -  $0(ftQ)  +  3kT  +  f(n0,T)  .  (229) 

However,  is  not  well  determined  so  we  calculate 

W  +  f(Q0’T)  “  E(Q0’T)  -  3kT  .  (230) 

The  resulting  points  are  plotted  versus  temperature  in  Figs.  46  a,  b, 
and  c.  We  fit  these  points  with  a  function  of  the  form 

VV  +  f(VT)  -  w  +  c<vt2  •  (231) 


with  results  shown  in  Table  IX. 


(fin,T)  (Ry) 


-0. 4500.  ♦50-0. 4554-0. 450-0, 450-0. 4551  -0.  OH 
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TABLE  IX 

COEFFICIENTS  IN  THE  FIT  TO  +  f  (fiQ,T) 

o  ,  3  W  C(V 

0U0;  _ (RjO _  _ (Ry/k2) _ 

232  -0.472520  +  1.3  x  10~4  (3.915  +  1.4)  x  10~9 

256  -0.473094  +  1.2  x  10~4  (2.598  +  0.75)  x  10~9 

270  -0.473215  +  0.9  x  l(f4  (3.630  +  0.72)  x  10-9 

The  values  of  C(£1q)  in  Table  IX  are  plotted  versus  volume  in  Fig. 

47*  We  fit  these  points  with  a  straight  line  and  arrived  at 

C(n0)  =  (0.657  +  11)  x  lo"8  +  (-0.126  +  4.3)  x  10~9  Hq  ,  (232) 

so  that,  in  the  same  manner  as  for  the  solid  bcc  and  hep  phases,  we  have 
specified  the  equation  of  state,  E(^,T),  for  the  glassy  and  liquid 
states  of  sodium  over  the  volume  and  temperature  ranges  of  our  calcula¬ 
tions  . 

C.  Melting 

We  here  discuss  in  more  detail  the  dynamics  of  melting  of  our 
molecular  dynamics  system.  We  concentrate  on  the  calculations  performed 

3 

at  an  atomic  volume  of  256  a^  shown  in  Fig.  45b. 

We  may  investigate  whether  or  not  a  particular  calculation  has 
melted  in  four  different  ways. 

First,  as  we  have  already  mentioned,  we  notice  which  curve  in  E,T 
space  it  lies — the  solid  or  the  liquid  curve.  The  point  labeled  1  in 
Fig.  45b  appears  to  be  on  the  solid  curve,  3  is  on  the  liquid  curve,  and 
2  is  somewhere  in  between. 

Second,  we  may  evaluate  the  atomic  distribution  of  particles  as  dis¬ 
cussed  In  Sec.  IV. D.  The  atomic  distributions  for  points  1,  2,  and  3  at 


£ 

1 


* 
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one  cime  are  plotted  as  the  solid  lines  in  Figs,  48a,  b,  and  c,  respec¬ 
tively.  From  these  plots  no  definitive  statement  may  be  made  about  the 
state  of  the  system.  However,  if  the  particle  positions  are  averaged 
over  time  for  point  1,  as  shown  by  the  dashed  curve  in  Fig.  48a,  the 
crystal  shell  peaks  are  seen  to  appear,  indicating  that  the  particles  are 
still  oscillating  about  their  lattice  positions.  When  the  same  averaging 
is  performed  for  point  3,  the  character  of  the  distribution  does  not 
change.  The  same  averaging  for  point  2  shows  less  peaking  than  1  but  more 
than  3.  Therefore,  we  have  another  definitive  difference  between  the 
liquid  and  solid  equation-of-state  points. 

Thirdly,  we  may  investigate  the  average  distance  that  a  particle  is 
from  its  starting  lattice  position.  We  do  this  by  calculating  the  mean- 
square  displacement,  Ar  ,  defined  as 

Ar2  -  <(r  -  r0)2>  ,  (233) 

where  r^  is  the  initial  position  of  the  particle  and  r  is  the  current 
position.  We  then  plot  Ar^  versus  time.  This  was  done  for  calculations 
of  points  1,  2,  and  3  in  Fig.  45b,  with  the  results  shown  in  Fig.  49. 

The  mean-square  displacement  for  point  1  has  settled  down  to  a  constant 
value  so  that  the  average  displacement  of  a  particle  from  its  original 
position  remains  constant.  The  calculated  value  for  point  3,  however, 
shows  a  reasonably  constant  slope,  indicating  that  atomic  diffusion  is 
indeed  occurring  with  particles  moving  away  from  their  original  lattice 
sites.  The  mean-square  displacement  for  point  2  shows  behavior  in  be¬ 
tween  these  two  extremes,  indicating  that  partial  melting  is  occurring. 

Finally,  we  may  look  at  plots  of  the  positions  of  the  particles 
in  space.  To  help  keep  track  of  particles,  we  have  chosen  to  connect 
nearest  neighbor  particles  within  each  close-packed  plane  by  lines. 
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Fig.  48. 

The  solid  lines  are  Che  atomic  distributions  at  one  Instant  and  the 
dashed  lines  are  the  distributions  of  the  particle  positions  aver¬ 
aged  over  200  cycles  for  points  (a)  1,  (b)  2,  and  (c)  3  of  Fig.  45b. 


The  initial  configurations  with  these  lines  drawn  have  already  been 
shown  in  Figs*  27  and  28.  When  a  particle  gets  more  than  8.0  a^  from  its 
near  neighbor,  the  lines  are  no  longer  drawn.  Figures  50a,  b,  and  c 
show  these  plots  for  points  1*  2,  and  3  of  Fig.  45b.  It  is  obvious  that, 
for  the  most  part,  the  particles  of  the  calculation  of  point  1  maintain 
their  positions  within  the  crystal  structure,  while  point  3  particles 
are  diffusing  through  the  system.  For  point  2  there  seem  to  be  regions 
where  diffusion  is  occurring  and  regions  where  it  is  not. 

From  the  four  techniques  mentioned  above,  we  can  state  that  point  3 
has  melted  and  point  1  has  not.  Point  2  is  in  some  intermediate  state. 

As  the  error  bars  on  this  point  in  Fig.  45b  indicate,  it  did  not  reach  a 
definite  equilibrium  state  after  a  rather  long  molecular  dynamics  calcu¬ 
lation.  The  time  history  of  the  temperature  of  this  calculation  is 
shown  in  Fig.  51.  We  notice  here  a  long  wavelength  oscillation  of  the 
average  temperature  is  present,  as  indicated  by  the  dashed  line.  This 
calculation  has  not  reached  equilibrium.  Point  3,  however,  has  reached 
equilibrium  on  the  liquid  curve  at  about  400  K  and  we  should  be  able  to 
identify  this  as  an  upper  limit  to  the  melting  temperature  of  sodium. 

We  may  make  no  such  statement  about  a  lower  limit  because  point  1  may  be 
metastable  in  the  solid  phase. 
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We  investigated  the  melting  of  sodium  at  atomic  volumes  of  232  a^ 

3  3 

and  270  a^  in  the  same  manner  as  was  done  for  the  256  a^  calculations 

described  above.  The  calculated  points  marked  with  arrows  in  Figs.  45a 
and  c  indicate  calculations  for  which  the  diffusion  and  thus  the  decrease 
in  temperature  was  proceeding  too  slowly  to  warrant  allowing  the  calcula¬ 
tions  to  proceed  to  equilibrium.  From  the  calculated  points  we  may. 


Fig.  50. 

Positions  of  the  particles  in  the  molecular  dynamics  system  viewed 
looking  down  the  x-axis.  Lines  have  been  drawn  between  near 
neighbors  within  a  close-packed  plane  at  the  initial  configura¬ 
tion.  a,  b,  and  c  show  the  positions  for  the  calculations  which 
produced  points  1,  2,  and  3  of  Fig.  45b,  respectively. 
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however,  place  upper  limits  on  the  melting  temperatures  as  430  K  for  the 

3  3 

232  a^  case  and  370  K  for  the  270  a^  case. 

We  compare  our  results  with  experimental  data  on  sodium  to  indicate 
how  realistic  our  molecular  dynamics  system  is.  The  estimated  upper 
limit  to  the  melting  temperatures  of  400  K  is  consistent  with  the  ob¬ 
served  melting  temperature  of  370  K. 

25 

Gingrich  and  Heaton  reported  the  experimental  atomic  distribution 
for  sodium  at  373  K.  It  is  shown  as  the  dashed  curve  of  Fig.  52,  where 
it  is  compared  with  our  calculated  curve  for  liquid  sodium  at  400  K. 

From  the  upper  curve  of  Fig.  49  we  calculate  the  self-dif fusion  co¬ 
efficient  for  liquid  sodium,  which  is  defined  as 


D 


(234) 


so  that  D  is  one-sixth  the  slope  of  the  curve.  We  calculate  this  slope 
2 

to  be  0.056  ag/T0,  which  yields  a  self-dif fusion  coefficient  of  D  * 

—5  2  26 

3.7  x  10  cm  /s.  The  experimental  value  reported  by  Faber  is  4.2  x 

10  ^  cm'Vs. 

We  also  compare  the  difference  between  the  solid  and  liquid  curves 
at  constant  temperature  with  the  latent  heat  of  fusion  for  sodium  which 

-3 

is  31.7  cal/g  -  2.32  x  10  Ry/ion.  The  difference  between  the  solid  and 

3  -3 

liquid  curves  at  400  K  for  the  256  a^  atomic  volume  case  is  1.7  x  10 
Ry/ ion. 

The  above  results  Indicate  that  the  molecular  dynamics  system  is 
reproducing  the  essential  characteristics  of  solid  and  liquid  sodium. 

D.  Dynamic  Phase  Change 

In  this  section  we  investigate  the  bcc-to-hcp  phase  change.  We 
have  already  mentioned  that  the  constraint  of  the  fixed  boundary 
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conditions  may  cause  a  system  to  be  metas table.  This  is  why  we  may  per¬ 
form  calculations  on  bcc  sodium  at  low  temperatures  even  though  the  hep 
system  is  the  preferred  phase.  Both  the  bcc  and  hep  672-particle  systems 
occupy  the  same  volume  but  the  shape  of  the  system  boxes  is  different,  as 
described  in  Sec.  IV. B.  At  the  end  of  Sec.  IV. B  we  mentioned  that  a  com¬ 
pression  in  the  y-direction  of  the  bcc  close-packed  plane,  to  form  a 
hexagonal  structure  and  a  corresponding  increase  in  the  x-  and  z-direc- 
tions  of  appropriate  magnitude  to  maintain  constant  volume,  may  allow  the 
hep  form  to  "fit"  in  the  correctly  shaped  box. 

We  made  this  change  of  shape  on  a  bcc  system  at  a  temperature  of 
50  BC.  The  calculation  was  allowed  to  equilibrate  for  150  time  units,  as 
shown  in  the  time  history  plots  of  the  potential,  kinetic,  and  total 
energies  of  Fig.  53.  The  atomic  distribution  of  the  system  at  this  time 
was  that  shown  in  Fig.  54a,  and  the  arrangement  of  particles  is  shown  in 
Fig.  55a.  For  this  figure  we  have  connected  the  nearest  neighbors  with 
lines  for  clarity.  (For  a  perfect  bcc  structure  without  thermal  motion, 
all  the  lines  would  be  either  horizontal  or  vertical  in  Fig.  55a). 

By  changing  the  shape  of  the  box,  we  created  a  body-centered  tetrag¬ 
onal  system,  thus  increasing  the  potential  energy  discontinuously  and 
doing  work  on  the  system.  The  jump  in  potential  and  total  energies  is 
evident  in  Fig.  53  and  the  body-centered  tetragonal  atomic  distribution 
and  particle  positions  are  shown  in  Figs.  54b  and  55b. 

The  system  was  allowed  to  equilibrate  on  its  own.  The  change  of 
shape  had  the  effect  of  pushing  the  system  over  a  potential  energy  bar¬ 
rier  and  allowing  the  particles  to  find  their  hep  configuration.  This 
is  attained  by  the  close-packed  planes  shifting  relative  to  one  another, 
which  is  characteristic  of  martensitic  phase  transitions.  This  shifting 
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Fig.  54. 

The  atomic  distributions  for  the  dynamic  phase  change  calculation 
for  (a)  the  bcc  system  before  the  change  of  shape  of  the  periodic 
box,  (b)  the  body-centered  tetragonal  unstable  system  immediately 
after  the  shape  change,  (c)  the  increased  temperature  hep  system 
after  the  spontaneous  phase  change,  and  (d)  the  cooled  hep  system. 


The  molecular  dynamics  system  viewed  down  the  y-direction  with  nearest 
neighbors  connected  by  lines  for  (a)  the  bcc  system  before  the  shape 
change,  (b)  the  body-centered  tetragonal  system  Immediately  after  the 
shape  change,  and  (c)  the  cooled  hep  structure  with  the  relative  plane 
positions  Indicated  by  A,  B,  or  C.  Outlining  the  figures  does  not 
represent  the  periodic  boundaries  and  Is  Included  for  reference  only. 
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is  indicated  by  the  arrows  in  Fig.  55b.  As  the  planes  shift,  the  poten¬ 
tial  energy  is  lowered  and  this  is  accompanied  by  a  corresponding  increase 
in  the  kinetic  energy,  thus  raising  the  system  temperature.  The  system 

reaches  equilibrium  as  shown  in  Fig.  53.  The  atomic  distribution  for 
this  elevated  temperature  system  is  shown  in  Fig.  54c.  We  then  cooled 
the  system  back  to  about  50  K.  The  atomic  distribution  of  Fig.  54d  shows 
the  peaks  characteristic  of  the  hep  structure.  The  particle  positions 
are  shown  in  Fig.  55c,  where  the  hexagonal  structure  is  apparent.  An  in¬ 
vestigation  of  the  stacking  of  the  planes  showed  that  the  system  did  not 
equilibrate  to  a  perfect  hep  structure,  which  is  an  ABABAB***  arrangement 
of  close-packed  planes  (as  discussed  in  Sec.  TV.B),  but  there  are  stacking 
faults  as  indicated  by  the  ABC  labeling  of  the  planes  of  Fig.  55c.  Such 
stacking  faults  are  prevalent  in  nature. 

We  performed  a  computer  experiment  on  a  similar  system  where  we 
heated  the  system  close  to  melting  and  cooled  it  down.  We  found  that  the 
system  would  return  to  the  hep  structure,  but  with  stacking  faults  in  dif¬ 
ferent  places. 

It  is  difficult  to  make  any  quantitative  statements  about  this  demon¬ 
stration.  However,  it  does  indicate  that  a  change  of  shape  of  the  peri¬ 
odic  box  is  necessary  to  observe  a  phase  change  of  the  type  demonstrated 
here.  A  shape  change  is  one  of  the  properties  that  characterizes  a  mar¬ 
tensitic  phase  change  in  nature.  A  change  of  shape  might  not  be  necessary 
if  the  system  size  were  very  large  or  if  free  volume  existed  which  would 
ease  the  constraining  effects  of  the  periodic  boundaries. 


159 


VII.  DISCUSSION  AND  CONCLUSIONS 

In  this  paper  we  have  demonstrated  that  the  molecular  dynamics 
technique,  coupled  with  an  interaction  potential  that  adequately 
describes  the  ion- ion  interaction,  can  be  used  to  study  the  macroscopic 
properties  of  a  simple  metal.  This  technique  is  unique  in  that  it 
provides  a  direct  calculation  of  the  anharmonic  terms  in  the  total  system 
energy.  We  calculated  and  presented  these  terms  as  a  function  of  volume 
and  temperature  for  the  solid  hep  and  bcc  phases  and  the  glassy  solid 
and  liquid  state  of  sodium. 

We  also  calculated  and  discussed  the  melting  of  sodium  and  were 
able  to  reproduce  the  experimental  heat  of  fusion  and  put  an  upper  limit 
on  the  melt  temperature.  We  showed,  by  comparison  with  the  experimental 
atomic  distribution  and  diffusion  coefficients,  that  our  calculations 
adequately  represent  the  dynamic  properties  of  sodium. 

As  a  final  demonstration  of  the  capabilities  of  this  technique, 
we  presented  the  results  of  a  calculation  that  reproduces  the  bcc-hcp 
martensitic  phase  transition  in  sodium  as  a  dynamic  process.  Although 
this  transition  was  artificially  induced  by  a  change  of  shape  of  the 
calculational  volume,  it  demonstrates  that  such  studies  are  feasible 
and  indicates  that  such  shape  changes,  which  occur  in  nature,  will  be 
a  necessary  part  of  future  studies. 

The  possibilities  of  future  work  using  this  technique  are  many. 
Although  we  have  the  ability  to  determine  the  system  energy  using  molec¬ 
ular  dynamics,  we  must  look  elsewhere  to  determine  the  system  entropy. 
This  may  be  done  by  performing  calculations  in  regions  where  the  theories 
are  known  to  work.  For  example,  a  quasi-harmonic  theory  may  be  used  to 


determine  the  entropies  at  low  temperatures  for  the  solid.  Ideal 
gas  or  hard-sphere  theories  may  be  used  at  high  temperatures  for  the 
liquid.  Once  the  entropy  is  known,  the  molecular  dynamics  results  may 
be  used  to  integrate  for  the  free  energies  and  thus  the  thermodynamic 
properties  will  be  completely  specified.  With  this  done,  the  phase 
change  and  melt  regions  may  be  determined  from  a  comparison  of  free 
energies . 

The  technique  should  also  prove  valuable  because  it  may  be 
extended  to  higher  density  and  temperature  regions.  Such  theoretical 
determination  of  the  equation-of-state  and  dynamic  properties  is  applic¬ 
able  to  many  areas  of  current  interest  such  as  the  study  of  shock-induced 
conditions.  The  fact  that  the  interaction  potential  is  volume  dependent 
means  that  regions  of  varying  densities,  which  exist  (for  example)  during 
the  shocking  of  a  material,  may  be  realistically  treated. 

However,  this  pseudopotential  theory  as  presented  here  is  limited  to 
simple  metals  and  compressions  of  about  50%  by  the  nearly-f ree-electron 
behavior  and  the  theoretical  constraint  that  the  ion  cores  do  not  overlap. 
The  extension  of  the  theory  to  handle  systems  with  more  complicated 
electronic  structures  is  being  done  (see,  for  example.  Refs.  1,  2,  6, 
and  7)  and  there  is  no  compelling  reason  why  the  theory  could  not  be 
developed  and  this  technique  applied  even  after  the  ion  cores  overlap 
and  ionization  occurs.  We  would  then  be  able  to  calculate  the  properties 
of  materials  at  very  high  temperatures  and  pressures  that  are  not  easily 
accessible  to  experiment. 
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APPENDIX  A 


CALCULATED  VALUES  OF 

THE  EFFECTIVE  ION- ION 

INTERACTION  POTENTIAL 

FORCE  FOR 

SODIUM  AT  AN  ATOMIC  VOLUME  OF  256  ag 

(The  potential  and  force  are  given  by  Eqs. 
plotted  in  Figs.  9  and  10.) 

(140)  and  (141)  and 

Radius 

Potential 

Force 

(ao) 

(Ry) 

(Ry/ag) 
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Radius 
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APPENDIX  B 


CALCULATION  OF  THE  KURTOSIS,  C,  FOR  A  MICROCANONICAL  ENSEMBLE 
OF  A  ONE-DIMENSIONAL  CHAIN  OF  HARD  RODS* 


For  a  microcanonical  ensemble  the  energy,  number  of  rods,  and  volume 
are  constant. 

E  -  energy 

P  -  radius  of  the  constant  energy  hyper sphere  in  momentum 
space  ■  (2mE) 

L  *  system  length 
N  *  number  of  rods 

*  momentum  of  ith  particle 

V  (p)  *  volume  of  an  N-dimensional  hypersphere  of  radius  p 
For  this  system  we  have 


E  •  T  ““  '  S  ■ 


N 

i-1 


u — 

P  *  ^  ^  p^  *  0  (center  of  mass  fixed) 
1-1 


pI1  “S  Pi  ’  n  *  l'2*' 

i-1 


Based  on  a  calculation  by  Brad  Lee  Volian,  Los  Alamos  National  Laboratory. 
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4 

$ 


J 

t* 


and 


UrnE)*5  <  P  <  [ 2m(E  +  AE)]1/2 


in  the  limit  as 


P  2  E  * 


29 


The  partition  function  for  the  microcanonical  ensemble  is 

ZN^N  /dPi  /dp2--'/dPN/dql/d<J2”-/d<1N  ’ 

where  *  implies  that  the  momentum  is  constrained  to  be  on  the  N-dimen- 

h 

sional  hypersphere  of  radius  P  *  (2mE)  ,  a  constant  energy  surface. 

L  -  NO  is  the  length  available  to  each  particle,  so  that 


N 

lim  [VM(P  +  AP)  -  VM(P)]  . 


N!h  AP 


The  average  of  the  quantity  P  is  given  by 


We  now  define 


r*  /-A 

rn(P)  -y  <*!-••/ dp Pin 

i-1 


so  that 


.  i. . 


p 


n 


11m 


I  (P  +  AP)  -  I  (P) 
n _ n 

V  (P  +  £f)  -  V  (P) 
n  n 
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I  (P)  is  an  N-dimensional  hypersphere  integral, 
n 

The  volume,  VN(P)  is  given  by 


VF) 


n 


N/2 


rt(N/2)  +  i] 


V1)P 


For  n  odd 

In(P)  *  0  , 


and  for  n  even 


in(p) 


nvnU> 


rf(N/2)  +  lirr(N  +  l)/2]  N+n 
r[(N  +  n)/2  +  l]r(l/2) 


We  now  expand 


(P  +  AP)N4tl  -  pn*>  -  ?^[(l  +  ff*  -  i] 

-  PN+Il[l  +  (N  +  n)  +  •  ‘  •  -  l] 


PN+n(N  +  n)  y-  +  ••• 


K 


so  that 


r  [  (N/2)  +  l)irr(N  +  l)/2]  .  _\„N+n  APj.  . 

In(P  +  AP)  -  In(P)  -  NVn(1)  F[(N  +  n)/2  +  X]r (1/2)  L  (N  n)  P 


and  similarly 


VN(P  +  AP)  -  VN(P)  -  VN(1)[(P  +AP)N  -  PN] 


NVn(1)PN^  + 


3 

1 


so  that 


*• 


v. 


p 


n 


11 m 


IN(P  +  £P)  -  IN(P> 
VN(P  +  SO  -  VN(P) 


,n  (N  +  n)  H  (N/  2)  +  1]T[(N  +  l)/2] 
T[(N  +  n)/2  +  l]r(l/2) 


and  we  see  that 


n  *  0 
n  -  2 
n  =*  4 


2  2 

=  P  ,  and 

7  ,  3PA . 


So,  on  using 


V 


n 


1  P 
N 


n 

n 


m 


and 


Pn 

Nmn 


(2mE)n/2  m  l/2E\n/2 

„_n  "  N'tn  / 

Nm 


1  /NkT\n/2 

N  V  m  / 


we  may  calculate  the  kurtosis,  C,  in  terms  of  the  velocity  moments 


kT 

m 


-4  _  3P4  _  „4 


+  2 


3N  / kT\2 
N  +  2  V  m  / 


so  that 


APPENDIX  C 


CALCULATED  EQUATION-OF-STATE  POINTS  FOR  hep,  bee,  AND  LIQUID  SODIUM 
The  equation-of-state  points  are  tabulated  according  to  the  ini¬ 
tial  configuration  of  the  system.  For  example,  the  high  temperature 
bee  calculations  are  reported  with  the  bcc  results  even  though  the  sys¬ 
tem  may  have  melted.  All  calculations  for  which  the  initial  condition 
was  the  glassy  state  are  reported  as  equation-of-state  points  for  liquid 
sodium.  The  reported  errors  are  the  standard  deviations  of  the  means 


and  represent  the  fluctuations  of  the  system  temperature  and  total  en- 


ergy  at  equilibrium. 

T 

Error 

Energy 

Error 

c 

(K) 

OQ 

(Ry) 

(x  10  °  Rv) 

hep  Structure 

98.83 

0.4 

-0.456570 

0.4 

=  232,  <t>  =  -0.458465 

198.01 

0.7 

-0.454657 

0.4 

(ERRATA:  Add  -0.015559  Ry 

to  energy  values) 

294.48 

0.5 

-0.452755 

0.6 

hep  Structure 

99.63 

0.3 

-0.456778 

0.2 

=  250,  =  -0.458682 

198.60 

0.5 

-0.454864 

0.3 

(ERRATA:  Add  -0.015960  Ry 

to  energy  values) 

293.98 

0.4 

-0.452962 

0.6 

hep  Structure 

10.03 

0.04 

-0.458394 

0.06 

=  256,  <J>  =  -0.458585 

49.68 

0.3 

-0.457634 

0.2 

(ERRATA:  Add  - 0.016089  Ry 

to  energy  values) 

99.19 

0.3 

-0.456685 

0.2 

198.11 

0.5 

-0.454771 

0.2 

293.27 

0.4 

-0.452873 

0.5 

V, 


hep  Structure 

=  270,  4>0  -  0.458131 

(ERRATA:  Add  -0.016374  Ry 
to  energy  values ) 


bcc  Structure 

-  232,  <J>0  -  -0.458447 

(ERRATA:  Add  -0.15559  Ry 
to  energy  values) 


bcc  Structure 

=  250,  $0  =  -0.458551 

(ERRATA:  Add  - 0.015960  Ry 
to  energy  values) 


bcc  Structure 

=  256,  *  -0.458447 

(ERRATA:  Add  -0.016089  Ry 
to  energy  values) 


170 

Error 


T 

(K) 

Error 

(K) 

Energy 

(Ry) 

(x  10'6  Ry) 

20.01 

0.08 

-0.457743 

0.1 

49.31 

0.3 

-0.457176 

0.2 

98.73 

0.3 

-0.456227 

0.2 

196.09 

0.4 

-0.454313 

0.3 

291.93 

0.4 

-0.452416 

0.3 

100.43 

0.5 

-0.456472 

0.4 

199.62 

0.8 

-0.454558 

0.6 

296.41 

0.8 

-0.452656 

0.6 

390.78 

0.6 

-0.450747 

1.0 

432.10 

0.5 

-0.447884 

0.6 

463.88 

-0.448818 

(not 

equilibrated) 

476.88 

-0.446919 

(not 

equilibrated) 

100.36 

0.6 

-0.456650 

0.4 

198.76 

0.8 

-0.454737 

0.5 

296.26 

0.4 

-0.452830 

0.6 

10.04 

0.05 

-0.458256 

0.07 

50.37 

0.2 

-0.457495 

0.1 

100.15 

0.6 

-0.456546 

0.4 

198.33 

0.8 

-0.454635 

0.5 

294.28 

0.3 

-0.452744 

0.3 

294.86 

0.3 

-0.452729 

0.4 

343.94 

0.4 

-0.451768 

0.4 

387.94 

0.6 

-0.450851 

0.5 

% 


399.13 


0.5 


-0.448906 


0.5 


T 

(K) 

Error 

(K) 

Energy 

(Ry) 

Error 

(x  10~6  Ry) 

bcc  Structure 

-  256,  $0  -  -0.458447 

408.32 

488.28 

0.6 

-0.449869 

-0.447009 

(not 

equilibrated ) 

0.7 

(cont) 

587.46 

1.0 

-0.445122 

1.0 

bcc  Structure 

=  270,  4>0  -  -0.457963 

99.11 

199.23 

0.5 

0.4 

-0.456066 

-0.454153 

0.4 

0.4 

(ERRATA:  Add  -0.16374  Ry 
to  energy  values) 

246.98 

0.4 

-0.453214 

0.5 

295.08 

0.3 

-0.452261 

0.5 

340.82 

0.5 

-0.451290 

0.6 

374.16 

-0.450357 

(not 

equilibrated) 

404.68 

-0.449395 

(not 

equilibrated) 

411.31 

0.7 

-0.448432 

1.0 

498.54 

0.6 

-0.446535 

0.7 

Liquid,  =  232 

10.01 

-0.456692 

(ERRATA:  Add  -0.015559  Ry 
to  energy  values) 

12.78 

19.24 

-0.456693 

-0.456538 

(not  su  le) 

31.79 

0.02 

-0.456269 

0.01 

82.54 

0.2 

-0.455298 

0.07 

125.17 

0.2 

-0.454516 

0.1 

183.37 

0.3 

-0.453373 

0.2 

248.42 

0.5 

-0.452092 

0.3 

329.61 

0.6 

-0.450440 

0.4 

403.75 


0.5 


-0.448639 


0.4 
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T 

00 

Liquid,  «  256  17.83 

(ERRATA:  Add  -0.016089  Ry  30.70 

to  energy  values) 

70.47 

137.61 

229.80 
229.96 
263.39 

358.78 

374.80 
482.73 

Liquid,  ^  =  270  7.10 

(ERRATA:  Add  -0. 016374  Ry  30.06 

to  energy  values) 

62.04 

142.79 
194.91 
262.31 
307.33 
396.01 


Error 

(K) 

Energy 

(Ry) 

Error 

(x  10-6  Ry) 

-0.456862 

(not  stable) 

0.03 

-0.456593 

0.03 

0.2 

-0.455727 

0.1 

0.4 

-0.454357 

0.1 

0.3 

-0.452499 

0.2 

0.4 

-0.452444 

0.3 

0.4 

-0.451672 

0.2 

0.3 

-0.449743 

0.2 

0.6 

-0.449574 

0.4 

0.9 

-0.447149 

0.7 

-0.456667 

(not  stable) 

0.04 

-0.456244 

0.02 

0.1 

-0.455619 

0.06 

0.4 

-0.454108 

0.2 

0.4 

-0.453033 

0.2 

0.6 

-0.451668 

0.3 

0.6 

-0.450628 

0.3 

0.4 

-0.448595 

0.3 

5 

■a 

vi 


>5 

,3 

4 

i 
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